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Part I 
Low Frequency Electromagnetic Modeling 
and Experiment: Tumor Treating Fields 


The Impact of Scalp’s Temperature A 
in the Predicted LMiPD in the Tumor giente 
During TTFields Treatment 

for Glioblastoma Multiforme 


Nichal Gentilal, Ariel Naveh, Tal Marciano, Zeev Bomzon, 
Yevgeniy Telepinsky, Yoram Wasserman, and Pedro Cavaleiro Miranda 


1 Introduction 


1.1 Tumor Treating Fields (TTFields) 


Tumor Treating Fields (TTFields) is an anti-mitotic cancer treatment technique 
used for solid tumors. It consists in applying an electric field (EF) with a frequency 
between 100 and 500 kHz in two perpendicular directions alternately to affect the 
mitosis of tumoral cells [1]. TTFields are FDA-approved for the treatment of recur- 
rent glioblastoma (GBM) since 2011, of newly diagnosed GBM cases since 2014, 
and of malignant pleural mesothelioma since 2019, after the promising results from 
the EF-11 [2], EF-14 [3] and STELLAR [4] clinical trials, respectively. The mecha- 
nisms by which TTFields act are not completely understood, but in-vitro studies 
showed that these EFs can slow down the rate at which tumoral cells divide or even 
completely arrest cell proliferation for EF intensities at the tumor of at least 1 V/ 
cm [1, 5]. 

In patients, TTFields are applied using a specific device named Optune 
(Novocure, Haifa, Israel) which consists in an electric field generator connected to 
four arrays with 9 transducers each. These arrays work in pairs to induce an EF in 
the tumor in two perpendicular directions alternately. The rationale behind the 
application of the fields in more than one direction is based on the results of the in- 
vitro study by Kirson et al. [5] that showed an increase of 20% in the number of 
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tumoral cells arrested during mitosis when two directions were used alternately 
compared to a single direction. 

Post-hoc analyses of clinical trials data showed that the time that a patient is on 
active treatment, i.e., the treatment compliance, significantly affects the expected 
overall survival (OS) [6]. Patients with recurrent GBM tumors that participated in 
the EF-11 trial and whose daily compliance was at least 18 hours had an OS of 
7.7 months, whereas those who were under treatment less than these 18 daily hours 
had an OS of 4.5 months. More recently, Ballo et al. [7] used data from newly diag- 
nosed GBM patients who took part on the EF-14 clinical trial to find a relation 
between TTFields dose and treatment outcome. The results showed that the OS 
survival is significantly different in patients whose average electric field in the 
tumor was at least 1.06 V/cm (24.3 months vs. 21.6 months), thus corroborating 
what was seen in in-vitro studies [1]. In the same study, the authors suggested a new 
metric, the power density (PD), that could be used to evaluate treatment efficacy 
instead of the electric field. By defining the dose in the tumor in the terms of PD, 
TTFields planning will be more like other cancer treatment techniques planning, 
such as radiotherapy, as this metric quantifies the energy that is deposited in tissues. 
The same study showed that the threshold that best divides patients into the two 
groups with the most significant difference in terms of OS is 1.15 mW/cm?. 


1.2 Optimization of the Array Placement: 
The NovoTAL System 


To optimize the PD in both directions, the array layouts are strategically placed on 
the scalp in regions that maximize the dose delivered to the tumor. The NovoTAL 
system (Novocure, Haifa, Israel) is an FDA-approved tool used to help certified 
physicians find the best array layouts for each patient. The importance of personal- 
ized array layouts is clearly shown in the literature by different authors. The study 
by Wenger et al. [8] was the first work to quantify the importance of creating per- 
sonalized treatment maps for each patient based on tumor location. In general, the 
average field intensity in the tumor increased between 32 and 45% when array posi- 
tioning was adjusted to the tumor location. Smaller improvements were seen in the 
work by Korshoej et al. [9], in which adapted layouts led to an enhancement of the 
average field strength in the tumor between 9% and 23%. 

In the NovoTAL system personalized head models are created and several 
array layouts are tested to find the one that yields the best option for each patient. 
One of the approaches that is most commonly considered comprehends an exten- 
sive search in which a finite number of layouts, available in the NovoTAL data- 
base, are tested and the one that induces the highest EF in the tumor is chosen. 
After the best layout is known small adjustments can be made to optimize the 
induced electric field. A detailed description of how this software works can be 
found in [10]. 
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1.3 Heat Transfer During Treatment 


In the studies mentioned above, the choice of which layout to use was based solely 
on the electric field delivered to the tumor. However, there are other factors that can 
affect the current that is injected in each array pair during the therapy and conse- 
quently the predicted treatment efficacy. As TTFields are applied for long periods of 
time, the temperature of the head increases inevitably because of the Joule effect. To 
avoid any thermal harm to tissues, the Optune device monitors the temperature of 
the scalp and keeps it around 39.5 °C by controlling the amount of current that is 
injected into each array pair. As the EF produced at the tumor by a given pair of 
arrays is proportional to the injected current, a decrease in the injected current 
implies a reduction in the EF at the tumor. 

Recent computational studies [11-13] modelled how TTFields were applied 
considering these thermal restrictions, quantified the temperature increases in tis- 
sues and predicted what physiological changes could occur when TTFields were 
applied. These works showed that heating was very localised, and it occurred mainly 
underneath the regions where the transducers were placed. The scalp was the tissue 
that heated up the most and the brain the one whose temperature varied the least. In 
each tissue, the temperature increased the most at the surface and quickly decreased 
with depth [13]. 


1.4 Aim of the Work 


Given that the scalp's temperature limits the dose that is administrated to the tumor, 
and consequently the predicted treatment effectiveness, in this work we investigated 
the impact of including this parameter on the choice of the best layout for TTFields 
treatment. 


2 Methods 


2.1 The Simplified Head Model (SHM) 


As a first step we used a simplified head model (SHM) to investigate the impact that 
the temperature might have in the metric used to choose the best layout. This model 
is composed by several spheroid shells each one representing a different tissue of 
interest (scalp, skull, CSF, brain and ventricles) as shown in Fig. 1. The dimensions 
assigned to each shell were based on the typical dimensions of the head of patients 
that participated in the EF-14 clinical trial. A spherical lesion with a radius of 
7.1 mm was also added to the model to mimic a GBM tumor. 
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Tissue a (mm) b (mm) 
Scalp 75 100 
Skull 68 93 
CSF 61 86 
Brain 59 84 

Ventricles 10 25 


4 
Fig. 1 Axial cut through the center of the model (plane z = 0). Each tissue was represented as a 
spheroid shell with the values of the main semi-axes, a and b, shown in the table on the right, in 
mm. The scalp is represented in orange, the skull in blue, CSF, which also fills the ventricles, in 


yellow and the brain in green. A spherical lesion, with the center at (34,18,0) mm and diameter of 
15.2 mm, in brown, was also added to the model to mimic a GBM tumor 


To maximize the number of layouts that could be tested we created several mod- 
els with transducer arrays placed in different regions of the scalp. All models were 
created using the Materialise 3-matic software. One important consideration to bear 
in mind when placing the arrays is that the 2 pairs should induce EFs with directions 
as close as possible to being perpendicular to increase the number of affected 
tumoral cells [1]. One pair, hereinafter referred to as pair 0, was placed as depicted 
in the top row of Fig. 2, in red. The array vector (AV), which was defined as the line 
that connects the two central transducers of each pair, passed through the center of 
the tumor. The direction of the array vector is 0° when it points along the y-axis and 
90° when it points along the x-axis. The perpendicular pair, pair 90, in blue in the 
center of the top row of Fig. 2, was rotated about the z-axis and its array vector also 
passed through the center of the tumor. In this work, we also considered the combi- 
nations between pair 0 and the pairs that were rotated by 75° and 105°, pair 75 and 
105, respectively, as it might not always be possible to ensure a complete perpen- 
dicularity between pairs in a real-case scenario. In total, we build seven different 
models in which the pairs were separated by 90° + 15° (Fig. 2). In all pairs the array 
vector passed through the center of the virtual lesion. 

The arrays are a 3x3 matrix of transducers separated by 44 mm in one direction 
and by 22 mm in the other, thus representing the Optune system. The radius of each 
transducer is 10 mm, and its thickness is 1 mm. Between each transducer and the 
scalp, a layer of gel with a radius of 10.3 mm and variable thickness was added to 
optimize the current flow into the head. 
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p —22 — 


Red: Pair O Red: Pair O Red: Pair O 
Blue: Pair 75 Blue: Pair 90 Blue: Pair 105 
Layout 0/75 Layout 0/90 Layout 0/105 


Y Y. A 


AAA 


Red: Pair 15 Red: Pair 30 Red: Pair 165 Red: Pair 165 
Blue: Pair 105 Blue: Pair 135 Blue: Pair 60 Blue: Pair 75 
Layout 15/105 Layout 30/135 Layout 165/60 Layout 165/75 


Fig. 2 Top view of the seven spheroid models used. For each pair there is a maximum of three 
possible complimentary pairs that can be used. For pair O these three pairs were rotated by 75° 
(upper row, left), 90° (upper row, middle) and 105° (upper row, right). For all layouts, the array 
vector of each pair passed through the center of the tumor, which is represented in black. The pair 
with the longest array vector is colored in red, whereas the pair with the shortest array vector is 
colored in light blue. The position of the tumor was the same for all layouts 


2.2 The Realistic Head Model (RHM) 


After the simulations with the simplified head model were completed, we used a 
realistic head model (RHM) to verify if the conclusions drawn hold for more com- 
plex geometries. The RHM was the same one used in previous studies on heat trans- 
fer during TTFields [11-13], This model was built based on the Colin27 template, 
and it was segmented into scalp, skull, CSF, grey matter (GM) and white matter 
(WM). Similar to what was done with the SHM, a spherical virtual lesion was added 
to mimic a GBM tumor. The latter consisted in a necrotic core, with a diameter of 
14 mm, surrounded by an active tumor shell, with a diameter of 20 mm. A descrip- 
tion of the complete framework followed to build this model can be found in [14, 
15]. The arrays were then placed on the scalp using Materialise 3-matic, in positions 
defined by the NovoTAL system. The final head model is shown in Fig. 3. 
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P A 
L 
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Fig. 3 Realistic head model used with the arrays placed in the regions suggested by the NovoTAL 
system. The anterior array (in blue) is paired with the posterior one (in green) and the left array (in 
red) is paired with the right one (in yellow). A: Anterior; P: Posterior; L: Left and R: Right 


2.3 Equations and Physical Properties 


As previously noted, current is injected alternately in two directions. Thus, to model 
TTFields, it is necessary to calculate the electric field distribution twice. At the fre- 
quency at which TTFields work for GBM treatment, 200 kHz, it is possible to use 
the electroquasistatic approximation of Maxwell's equation. In this case, the calcu- 
lation of the electric potential is simplified, and Laplace's equation can be used: 


V[(o+ je,2x f) V6 ] - 0 (1) 


where o is the electric conductivity (S/m), j the imaginary constant, e, the absolute 
permittivity (F/m), f the frequency (Hz) and ¢ the electrostatic scalar potential (V). 

As Optune controls the potential difference between the two arrays of the same 
pair, we imposed a Dirichlet boundary condition (V 2 + Vo) at the outer surfaces of 
the transducers of the same array. At the remaining outer boundaries of the model, 
we assumed a null value for the normal component of the current density, whereas 
at the internal boundaries we considered continuity of the latter quantity. The elec- 
tric field distribution was obtained for the two array pairs in each layout. 

After solving the electric studies, we used Pennes' equation to solve for the tem- 
perature distribution: 


oT 


pe =V{KVT) + ope (T, -T) + Qn +0 lel Q) 


met. 
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where p is the density (kg/m), c the specific heat at constant pressure (J/(kg °C)), T 
the temperature (°C), t the time (s), k the thermal conductivity (W/(m °C)), @ the 
blood perfusion (1/s), Qme: the metabolic rate (W/m?), assumed to be constant in the 
simulations, and Æ the electric field vector (V/m). The last term on the right-hand 
side represents the contribution of the EFs in increasing the temperature. Its values 
were taken alternately from each array pair. The subscript b stands for blood and T, 
was 36.7 °C in all studies. 

At the outer boundaries, we accounted for energy exchanges with the environ- 
ment, at 24 °C, through convection and radiation, whereas at the internal boundaries 
we assumed that heat conduction occurred between adjacent tissues and materials. 

The values assigned to each physical property were based on the literature. The 
electric parameters were retrieved from Ballo et al. [7] and the thermal parameters 
were the same as the standard values reported in Gentilal and Miranda [12]. 

All simulations were performed in COMSOL Multiphysics v5.2a. 


2.4 Conditions of the Simulations and Metrics Used 


As described in the work by Chaudry et al. [16], during treatment planning with 
NovoTAL the electric field in the tumor is predicted by injecting 900 mA of current 
(amplitude) into each array pair independently. Additionally, in the work by Ballo 
et al. [7], the local minimum power density (LMiPD) is defined as the lowest of the 
two power densities induced in each voxel in the tumor. Thus, we combined these 
two approaches and, when only the electric field was considered as a metric to 
evaluate the effectiveness of a layout, we injected 900 mA into each array pair, and 
we calculated the LMiPD in the tumor. As noted by Wenger et al. [10], it is desirable 
to use MR images with a resolution of at least | mm x 1 mm x 1 mm for TTFields 
planning. Thus, the values of the power density in the tumor were exported from 
COMSOL considering a regular grid spaced by 1 mm in every direction. After the 
minimum power density was obtained for each voxel, we calculated the average 
LMiPD in the tumor for each layout. Each electric study took around 4 minutes 
using the SHM and 30 minutes using the RHM in a workstation with dual core Intel 
Core 19-10900X X-series processors clocked at 3.7 GHz and with 64 GB of RAM. 

When the temperature was also accounted for during treatment planning a differ- 
ent approach had to be followed as the amount of current that could be injected into 
each array pair depended on the temperature of the scalp. We used a method previ- 
ously reported [13] in which we started by injecting 400 mA amplitude of current 
into each pair alternately with a switching time of one second. We then predicted the 
maximum temperature that the scalp would reach underneath each array pair in 
regions where the thermistors are typically placed, and we iteratively fine-tuned the 
injected current to induce the desired maximum temperature underneath both pairs. 

Predicting the temporal evolution of the temperature distribution required a time- 
dependent study, with the electric field in Eq. (2) switching between that produced 
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by each of the two array pairs every second. Only the first 5 minutes of treatment 
were simulated due to the time taken by these simulations. In order to estimate the 
steady state temperature, we fitted the following general equation to the variation of 
the highest temperature in the scalp: 


T - C, (1- exp(-t/C;)) C, (3) 


where C, (°C) represents the maximum contribution that the EF has in increasing 
the temperature of the scalp, C; (s) is a time constant related to how quick the heat 
is transferred and C; (°C) represents the initial temperature of the scalp. 

The maximum temperature, 7"^', can then be predicted by taking the limit as t 
tends to infinity: 


T™ = lim(r —o)T =C, +C, (4) 


All data curve fitting were performed using Matlab’s 2020 curve fitting toolbox 
(cftool). 

In the case of the RHM the desired temperature limit was the optimal working 
point for the Optune device, 39.5 °C. For the SHM, we had to slightly increase this 
threshold, to 40.5 ?C, as the minimum current that is typically injected in patients, 
400 mA, already induced a temperature higher than 39.5 °C. The time-dependent 
studies took around 10 hours for the SHM and 34 hours using the RHM in the work- 
station aforementioned. 


3 Results and Discussion 


3.1 Simplified Head Model 
3.1.1 LMiPD Values When Only the Electric Field Is Considered 


The values of the LMiPD when 900 mA (amplitude) were injected into each pair 
using the simplified head model are presented in Table 1. 

The values of the LMiPD were very similar for all layouts (range: 6.30-6.88 mW/ 
cm?). The worst layout for this tumor position, layout 165/75, induced a LMiPD that 
was only 896 lower than the one induced by the best layout, layout 15/105. These 
values were mainly limited by the lower EF induced by the pair with the longest 
array vector. In general, a longer array vector corresponds to a more resistive current 
path. However, resistance values also did not change considerably across layouts: 
90-93 Q for the pair with the highest AV and 77-82 Q for the complimentary pair. 
Under these conditions, the parameter that defined which was the best layout is the 
distance between the tumor and each array (Fig. 2). In layout 15/105 the tumor is 
very close to two arrays of different pairs thereby increasing the induced PD in both 
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Table 1 LMiPD values in the tumor when 900 mA (amplitude) were injected into each pair. Each 
layout is named after the two pairs which compose it, where the first pair is the one with the longest 
array vector and consequently the highest resistance 


Injected current (mA) Injected current (mA) LMiPD 
Layout First pair Second pair (mW/cm?) 
0/75 900 900 6.63 
0/90 900 900 6.61 
0/105 900 900 6.57 
15/105 900 900 6.88 
30/135 900 900 6.71 
165/60 900 900 6.31 
165/75 | 900 900 | 6.30 


directions. The distance between the tumor and the closest array of pair 15 is 55 mm, 
a value only slightly higher when compared to the distance to the closest array of 
pair 105, 48 mm. On the other hand, in layout 165/75, the distance between the 
tumor and the closest array of pair 75 is only 36 mm, but this value increases to 
78 mm when the distance to the closest array of pair 165 is considered. Thus, the 
LMiPD in the tumor will be higher for layout 15/105 compared to layout 165/75 for 
the same injected current. 

Based only on the LMiPD values, the layout that would be chosen for this model 
is layout 15/105. 


3.1.2 LMiPD Values When the Temperature Is Accounted for 


To investigate how the temperature impacts treatment planning we started by inject- 
ing 400 mA into each array pair alternately with a switching time of 1 second and, 
at the end of the simulation, we predicted the maximum temperature that the scalp 
would reach underneath both pairs. We then iteratively fine-tuned the injected cur- 
rent to induce the highest average current, considering both pairs, that would lead to 
a maximum of 40.5 °C in the scalp. 

Taking layout 0/75 as an example, Fig. 4 shows the maximum temperature varia- 
tion of the scalp for different sets of injected currents. 

The maximum temperature in the scalp for each set of injected currents was pre- 
dicted using Eq. (4) and it is presented in Table 2. The minimum adjusted-R? (A-R?2) 
obtained during the curve fitting process was 0.9995. 

An increase of 200 mA of current in pair 75, from 400 mA to 600 mA, led to a 
temperature increase of 2.8 °C in the scalp underneath this pair and of 0.4 °C under- 
neath pair O as the two pairs were placed very close to each other. A more detailed 
description of the impact of the distance between array pairs is provided in the next 
subsection. 

The same rationale was followed for the remaining 6 layouts. The LMiPD was 
then quantified for the current values that led to 40.5 ?C in the scalp (Table 3). The 
lowest A-R? was 0.9994. 
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Layout 0/75: underneath pair 0° 35 Layout 0/75: underneath pair 75° 


Maximum temperature (°C) 
Maximum temperature (°C) 


0 50 100 150 200 250 300 
Time (s) Time (s) 


Fig. 4 Scalp's maximum temperature variation underneath pair O (left) and pair 75 (right) for dif- 
ferent values of injected current. Each line corresponds to a different set of values denoted by X/Y 
mA, where X is the current injected in pair 0 and Y the current injected in pair 75. Current was 
injected alternately into each pair with a switching time of 1 second. Changing how much current 
was injected in pair 75 also changed the temperature of the scalp underneath pair 0 


Table 2 Maximum temperature predicted in the scalp underneath each pair of arrays when 
different sets of current were injected in each pair alternately with a switching time of 1 second 


Injected current (mA) Maximum scalp's temperature Maximum scalp's temperature 
Pair 0/Pair 75 Underneath pair 0 (°C) Underneath pair 75 (°C) 
400/400 40.1 37.7 

400/550 40.4 39.5 

400/600 40.5 40.5 


Table3 LMiPD values in the tumor for each layout when current was injected into each array pair 
alternately with a switching time of 1 second. In every case, the amount of current injected was 
iteratively refined to induce a maximum temperature of 40.5 °C in the scalp underneath each pair 


Injected current (mA) Injected current (mA) T + LMiPD 
Layout First pair Second pair (mW/cm?) 
0/75 400 600 1.31 
0/90 410 625 1.38 
0/105 423 600 1.46 
15/105 400 575 1.44 
30/135 425 600 1:71 
165/60 423 600 1:39 
165/75 423 625 1.39 


Accounting for the temperature led to a variation of 31% in the LMiPD values 
(range: 1.31-1.71 mW/cm?), a value higher than the 8% reported when only the 
electric field was used as a metric to predict treatment effectiveness. In general, the 
current injected into the most resistive pair did not vary by more than 25 mA (range: 
400-425 mA), whereas the current injected into the complimentary pair varied 
between 575 and 625 mA. 
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Based on Table 3, the best treatment option is layout 30/135. When only the 
electric field was used as a criterion to choose layouts (Table 1), this layout would 
be the second choice, after layout 15/105. On the other hand, layout 15/105 would 
now be the third choice if information regarding scalp’s temperature was included 
during treatment planning. 


3.1.3 Comparison Between LMiPD and T + LMiPD 


Figure 5 depicts the variation in the ranking of the layouts based on the values of the 
LMiPD when the temperature is not accounted for during treatment planning 
(LMiPD) vs. when it is (T + LMiPD). 

Based on these data, it is noteworthy that layouts 0/75, 0/90 and 15/105 dropped 
in rank when information regarding scalp temperature was included in treatment 
planning. By definition, the LMiPD is limited by the pair that induces the lowest 
power density in the tumor. As less current is injected into the most resistive pair in 
these three layouts compared to the other four (first pair, Table 3), the LMiPD will 
be reduced. As it can be seen in Fig. 2, some pairs in these three layouts are close to 
each other and thus the current that is injected into one pair will not only increase 
the temperature in scalp regions underneath that pair but also contribute to increas- 
ing the temperature underneath the complimentary pair. This can be seen quanti- 
tively in the data presented in Table 2, where an increase of 200 mA in pair 75 led 
to an increase of 0.4 °C underneath pair 0. 

In particular, two pairs in layouts 0/75, 15/105 and 0/90 created a temperature 
hotspot near the two closest transducers of different arrays, as shown in Fig. 6. In 
the first two layouts the shortest distance between transducers of two different pairs 


—9—0/75 —9— 0/90 —— 0/105 —*— 15/105 — — 30/135 —®— 165/60 —9— 165/75 


1 I mdi as ER ITA 


2 


best; 7=worst) 


5 Mun — ca"! 


Ranking order (1 


LMiPD T*LMiPD 
Criterion 


Fig. 5 Variation in the ranking of the layouts based on two different criteria: only the electric field 
in the tumor (LMiPD) or the electric field in the tumor plus information regarding the temperature 
of the scalp (T + LMiPD). Each line represents a different layout 
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Fig. 6 Temperature distribution at the end of the simulation (t = 5 minutes) in layouts 0/75 (left), 
15/105 (middle) and 0/90 (right), when 400/600 mA were injected in the first, 400/575 mA in the 
second and 410/625 mA in the third. Due to the proximity of some transducers of different pairs in 
these three layouts temperature hotspots occurred (black arrows) which limited how much current 
could be injected into each pair. In the first two the shortest distance between transducers of differ- 
ent pairs was 2.2 mm, whereas in the third it was 10.0 mm. Temperature scale is in °C 


was 2.2 mm. In layout 0/90, this distance was 10.0 mm, which created a weaker but 
still significant hotspot. 

This highlights the importance of considering scalp's temperature during 
TTFields planning, especially in patients with smaller heads in which there is a high 
probability for the transducers of different pairs being very close to each other. 

The average value of the LMiPD across layouts decreased from 6.59 mW/cm’, 
when only the electric field in the tumor was considered, to 1.51 mW/cm?, when 
information concerning the temperature of the scalp was also included, a reduction 
of 77% on the average LMiPD predicted in the tumor. This is due to the large reduc- 
tion in injected current imposed by the limit on the scalp temperature. 

This simplified head model allowed us to gain important insights on some of the 
most important parameters when planning TTFields. However, the conclusions 
drawn from it might only be useful qualitatively due to some limitations such as the 
small variation in the values of the metrics used and the fact that the temperature 
threshold had to be increased from 39.5 °C to 40.5 °C. The latter might suggest that 
the heat mechanisms, mainly the Joule heating term, is overcontributing to the tem- 
perature increases as this threshold was reached for low values of injected current 
compared to what is typically seen in patients. In the next section, we present a 
preliminary study where we performed the same type of analysis using a realistic 
head model. 


3.2 Realistic Head Model 
3.2.1 LMiPD vs. T + LMiPD 
In the realistic head model, the LMiPD in the tumor was also first calculated by 


injecting 900 mA into each pair, which yielded 0.81 mW/cm?. This value is below 
the 1.15 mW/cm*, reported by Ballo et al. [7], which is the threshold that best 
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divides GBM patients into two groups with the most significant differences in the 
overall survival. One important factor limiting the LMiPD in this head model is the 
electric field produced by the anterior-posterior (AP) configuration. The power den- 
sity in the tumor induced by this pair was 0.88 mW/cm’, a value far below the 
3.32 mW/cm? induced by the left-right (LR) configuration. The virtual lesion was 
placed in the right hemisphere of the brain, close to the lateral ventricles but mid- 
way between the anterior and posterior arrays, where the electric field produced by 
the AP pair was much lower than to the one induced by the LR pair. In general, 
tumors located in deeper regions of the brain, as in this case, require a more robust 
placement of the layouts as clearly shown in the work by Korshoej et al. [9]. 
According to the results of the latter work, a layout that might yield a higher LMiPD 
in our head model could be one in which all pairs are maintained at the same height, 
but rotated by 45°. This would decrease the power density induced by the LR pair 
but increase the one produced by the AP pair, which is the pair that is strongly limit- 
ing the LMiPD. 

We then included information regarding the temperature of the scalp in the simu- 
lations. We followed the same approach as before and our results indicated that if 
580 mA were injected into the AP pair and 850 mA in the LR pair alternately, the 
maximum temperature that the scalp would reach underneath both pairs would be 
39.5 °C. The lowest A-R? obtained for the data used to drawn these conclusions 
was 0.998. 

More current could be injected in the left-right configuration because the head is 
less resistive in this direction than in the AP direction. Figure 7, shows the tempera- 
ture distribution in scalp’s surface at the end of the simulation (t = 5 minutes). As 
expected, heating was very localised, and it occurred mainly underneath the regions 
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Fig.7 Scalp surface temperature after the first 5 minutes of treatment when 580 mA of current 
were injected into the AP pair and 850 mA into the LR pair, alternately with a switching time of 
1 second. The maximum temperature predicted in the scalp under these conditions is 39.5 °C 
underneath both pairs. As it can be seen, heating is very localised, and the temperature increases 
underneath one pair do not significantly affect the increases underneath the complimentary pair. 
Temperature scale is in *C 
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where the arrays were placed. Contrary to what occurred with layouts 0/75, 15/105 
and 0/90 in the SHM (Fig. 6), in this model the two pairs were relatively far from 
each other, and thus common temperature hotspots did not occur. 

The LMiPD was then calculated for these conditions (T + LMiPD) and com- 
pared with the values obtained when only the electric field was used as a criterion. 
The LMiPD induced by the AP configuration decreased from 0.88 to 0.36 mW/cm? 
and the one produced by the LR pair decreased from 3.32 to 2.96 mW/cm?. The 
overall LMiPD predicted in the tumor decreased from 0.81 mW/cn? to 0.36 mW/ 
cm?, when the temperature was also included in treatment planning, a 5696 reduc- 
tion. This suggests that taking the scalp's temperature into account is likely to intro- 
duce large variations in the LMiPD also when using realistic head models. 


4 Conclusions and Future Studies 


The results obtained using the SHM model suggested that the best treatment layout 
in terms of LMiPD was not necessarily the best one when the temperature was 
included in treatment planning. In some cases, layouts in which two array pairs 
were close to each other yielded the highest average EF in the tumor, but the occur- 
rence of a common temperature hotspot significantly limited the current that could 
be injected into each array pair, and thus the predicted treatment efficacy when the 
temperature was accounted for. However, adding this additional step during 
TTFields treatment planning might be very time-consuming. In our workstation, it 
took 10 hours to compute each time-dependent simulation. The process of fine- 
tuning the values of the injected current and of testing several layouts for each 
model might lead to computational times that are unfeasible from a practical point 
of view. Thus, further investigation is needed to include thermal studies in the treat- 
ment planning workflow. One important consideration that might be helpful is to 
define a minimum distance between the arrays to prevent the occurrence of tem- 
perature hotspots. 

The results obtained with the RHM are a part of a preliminary study that is under 
investigation. The data presented here were intended to show that the LMiPD pre- 
dicted in the tumor can also significantly decrease in realistic head models when the 
temperature of the scalp is considered during treatment. We are currently evaluating 
the impact of these decreases using different array layouts obtained through the 
NovoTAL system. In some of them, transducers of different pairs are also very close 
to each other, as it occurred with the SHM, and thus significant decreases in the 
predicted LMiPD are expected. Furthermore, we are also investigating different 
metrics that could be used alongside LMiPD to predict more realistically treatment 
effectiveness. 
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1 Introduction 


Tumor Treating Fields is a well-established treatment modality for glioblastoma 
(GBM) [1]. One factor that influences the efficacy of TTFields therapy is the elec- 
tric field strength. High electric field intensities reduce the rate of tumor cell divi- 
sion in vitro [2] and increase progression-free survival (PFS) and overall survival 
(OS) in GBM patients, thus demonstrating the TTFields therapy dose-response rela- 
tionship [3]. 
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Korshgj et al. proposed in 2016 a novel surgical technique “skull-remodeling 
surgery” (SR-surgery) to increase the TTFields intensity. SR-surgery involves 
removing bone from the skull in the form of craniectomy, burrholes, or skull thin- 
ning. The hypothesis is by removing parts of the cranium, the resistance created by 
the bone is reduced and the electric field intensity is strengthened focally in the 
underlying tissue. By using realistic computational head models and finite element 
method it was concluded that SR-surgery could potentially increase the TTFields 
intensity by 60-70% in superficial tumors. Furthermore, it was concluded that sev- 
eral small burrholes would induce a greater increase in field strength than one large 
craniectomy per skull defect area, which is an important clinical safety consider- 
ation [4]. Computational studies also indicated that TTFields electrode array place- 
ment could be optimized when taking SR-surgery into account [5]. 

Feasibility of the SR-surgery concept was recently demonstrated in a small 
investigator-initiated, single-center, open-label phase 1 trial (OptimalTTF-1) from 
2016 to 2019 with 15 first recurrent glioblastoma patients. For ethical- and risk-to- 
reward reasons one of the inclusion criteria was that the SR-surgery should increase 
the TTFields intensity by a minimum of 25%. This was validated using computa- 
tional methods on patient-specific models. The median relative increase in the field 
intensity was 32% (range 25;59). The trial concluded that SR-surgery was safe 
without additional toxicity [6]. 

Subsequently, a randomized, comparative, multi-center, investigator-initiated 
interventional phase 2 trial (OptimalTTF-2) was launched in October 2020 [7]. The 
trial will include 84 first recurrence GBM patients randomized 1:1 to TTFields ther- 
apy with or without SR-surgery. All patients will receive physician’s choice medical 
oncological therapy. The primary endpoint is overall survival. 

To simplify and standardize the intervention, the SR-surgery procedure is con- 
ducted using a standardizing operating procedure (SOP) that is easy to understand 
and replicate with minimum required training. 

The specifications of the SOP are given below, and the following work shows 
preliminary documentation for the approach. 


2 Method 
2.1 Standardizing the SR-Surgery Configuration 


A five burrhole SR-surgery configuration was chosen as shown in Fig. 1, based on 
previous research indicating higher efficacy compared to complete craniectomy [4]. 
The general rationale behind this configuration was to induce significant field 
enhancement in the peritumor and tumor region without compromising patient 
safety and with no requirement for protective headgear. 
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Fig. 1 Five burrholes each 15 mm in diameter. Central burrhole placed above the tumor 
(marked in red) 


Suggested clinical procedure: 


1. Five burrholes are created each 15 mm in size. 

2. The central burrhole should be placed directly above the tumor resection cavity 
or tumor remanent. 

3. Electrodes are placed so that edge electrodes from one array in a pair are located 
above the holes or very close to the holes. The other array in the pair should be 
placed directly opposite on the other side of the head. Preferably, a direct line 
between the transducer located above the burr-holes and the central transducer 
on the opposite array in the pair should pass through the region of interest, i.e., 
the residual tumor or a region of the resection cavity with a predicted high risk 
of recurrence. This configuration should be achieved from both TTFields pairs. 


2.2 Technical Aspects of the SR-Surgery 


To support the neurosurgeons in implementing the procedure outline above and to 
ensure uniformity across the phase II trial centers, 3D printed, and sterilized tem- 
plates were made to guide the neurosurgeons during surgery. Four templates with 
different radius curvatures (0.0, 3.2, 6.4, and 9.6 mm) were included in each steril- 
ized pack to improve fitting in different areas of the skull (Figs. 2 and 3). 

The neurosurgeon will decide perioperatively, after completion of the tumor 
removal where the bulk of the resection cavity or tumor remanent lies and place the 
central burrhole directly on top. If no residual tumor is left, the surgeon would place 
the central burrhole above the resection cavity or a part of the surrounding brain 
tissue where future recurrence is expectedly likely to occur. Placement of the bur- 
rholes does not need to be restricted to the craniotomy existing bone plate used to 
access the tumor during resection. The template is placed as described above and 
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Fig. 2 The four 3D 
printed templates are 
shown in a sterilized 
peel-pack. Each template 
curves slightly more than 
the next 


Fig. 3 The concept and schematic of a 3D printed template 


the five burrholes are marked with a high-speed electrical drill. The template is then 
removed during drilling and repeatedly repositioned until the right configuration is 
achieved. The additional procedure time is approximately 5-10 min, with minimal 
risk of dural- and cortical damage. Two examples are given in Figs. 4 and 5. 


2.3 Validating the SR-Surgery Configuration via Simulations 


To validate this SR-surgery configuration we used simulations of the electric fields 
generated by the TTFields therapy using a detailed head model constructed from 
structural MR images and the Finite-Element Method (FEM). The head model was 
initially created from a dataset of a healthy participant which was then adapted to 
emulate a trial patient's pathology based on their postoperative MRI and CT of 
the head. 
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Fig. 4 A post-operative 
CT of the skull in a patient 
that underwent SR-surgery. 
In this case, the bone flap 
was big enough to fit all 
five burrholes while 
simultaneously having the 
central burrhole right 
above the resection cavity 


Fig. 5 Perioperative 
photo. It shows a bone flap 
that is big enough to fit all 
five burrholes, however, 
the entire configuration 
was moved so that the 
central burrhole was over 
the bulk of the resection 
cavity. Therefore a “half 
burrhole” was created on 
the bone flap and the other 
half on the skull 


More specifically, we created the computational head model using the example 
dataset “Ernie” in the SimNIBS software package as a proof-of-concept demonstra- 
tion. The “Ernie” dataset corresponds to a young and healthy subject, including a 
high-resolution T1 and a T2 weighted image. Details on the “Ernie” dataset can be 
found in the documentation from the SimNIBS toolbox. 

In our first step, an automated tissue segmentation was performed on the T1 and 
T2 weighted images of the "Ernie" dataset. The initial segmentation includes the 
following eight tissue compartments: white matter, gray matter, CSF, scalp, skull, 
muscle, blood, and eyeballs. 
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In our second step, the segmented image was visually inspected, and the configu- 
rations of the virtual tumor resection cavity, the residual tumor, and the burrholes 
were set manually based on the trial patient’s postoperative MRI and CT of the head 
and were as follows: 


1. A virtual sphere-shaped tumor resection cavity with a 2.5 cm diameter in the 
“Ernie” head model (Fig. 6a). 

2. A cylinder-shaped funnel on the top of the tumor resection cavity mimics the 
surgical entry to the tumor. The funnel track has a 0.8 cm diameter. 

3. A sphere-shaped tumor remanent of a 2.5 cm diameter underneath the tumor 
resection cavity. 

4. Virtual SR-surgery was applied to the head model. The SR surgery entailed plac- 
ing five burrholes of 1.5 cm diameter on the skull with the central burrhole 
directly above the resection cavity. We investigated five different configurations 
of burrholes as described in Sect. 2.5. 


In our third step, the head models were created using the segmented image with 
virtual tumor cavity, residual tumor, and with or without burrholes. The final head 
mesh with each burrhole configuration consisted of a total number of approximately 
4,750,000 tetrahedral elements, assigned to more than 11 types, including while 
matter, gray matter, CSF, scalp, skull, muscle, blood, eyeballs, tumor resection cav- 
ity, residual tumor, and burrholes. 

In our fourth step, the electrode arrays were placed on the skin surface for each 
head model (see Fig. 7) and conductivity values were assigned to the different com- 
partments, consisting of skin (0.25 S/m), bone (0.010 S/m), CSF (1.654 S/m), grey 
matter (GM) (0.276 S/m), white matter (WM) (0.126 S/m), residual tumor 0.24 S/m) 


(b) 


Fig. 6 (a) The head model is based on a healthy person “Ernie” MRI and the following spherical 
2.5 cm resection cavity (blue) and spherical 2.5 cm residual tumor (yellow) was manually added 
based on a trial patient’s postoperative MRI and CT. (b) shows the five burrholes, each 1.5 cm in 
diameter. The central burrhole is placed directly over the resection cavity 
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(a) 0° (b) 45° (c)90° = (d) 135°  (e)180? 


Fig. 7 Surface reconstruction of the head model showing different electrode array positions used 
for simulation. Electrode arrays are comprised of 9 electrodes of 20 mm in diameter arranged in a 
3x3 array pattern. The center-to-center distances between the electrodes are 45x22 mm. A 
15-degree stepwise rotation of an electrode array was conducted around a central craniocaudal axis 
in the same horizontal plane, corresponding to degrees of 0, 15, 30, 45, 60, 75, 90, 105, 120, 135, 
150, 165, and 180 for a total of 13 different positions. Figure (a)-(e) give five exemplary positions. 
Electrodes are paired gray with white 


and necrotic tissue (1.0 S/m). We created a custom version of SimNIBS and pro- 
vided scripts for automated simulation of TTFields induced electric fields for varia- 
tions of the electrode array positions (described in the following section). 

At last, the results of the simulations were visualized using Gmsh. 


2.4 Electrode Array Placement 


Previous research indicates that SR-surgery alters the distribution of the field den- 
sity in the skull and that there are electrode array positions that are the most optimal 
when taking into consideration the tumor- and SR-surgery positions [5, 8]. 

To systematically examine the optimal placement of the electrode arrays a 
15-degree stepwise rotation of an electrode array was conducted around a central 
craniocaudal axis in the same horizontal plane, corresponding to degrees of 0, 15, 
30, 45, 60, 75, 90, 105, 120, 135, 150, 165, and 180 for a total of 13 different posi- 
tions as shown in Fig. 7. 

FEM calculations were performed for each electrode placement with and with- 
out the chosen SR-surgery configuration. 


2.5 SR-Surgery Location in Relation to the Tumor 
on the Electric Field 


To test the hypothesis that the burrholes should be placed directly above the tumor 
resection cavity to maximize the electric field directly underneath, three additional 
simulations were performed. Each SR-surgery configuration was moved to three 
alternative positions as shown in Fig. 8. The electrode arrays were rotated as 
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wu 


(d) Burrholes moved *up and back" 


(e) Without burrholes 


Fig. 8 (a) Top left corner shows the original placement directly above the resection cavity. The 
three locations were moved 3 cm in each direction “up” (b), “back” (c), and “back and up" (d). (e) 
is with intact skull that serves as a control 


described in Fig. 7 while the electric field was calculated. For all SR-surgery, posi- 
tions electrodes array rotation overlapped the burrholes. These simulations were 
performed to examine if small variations in the SR-surgery location had an impact 
on the electric field in the resection cavity and residual tumor. 
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3 Results 


1. Burrholes significantly increased the mean, median, and peak electrical field 
intensities in residual tumor tissue and the resection cavity (Figs. 9 and 11). 

2. The effects of burrholes on the distribution of field strengths in healthy tissues 
were minor (Figs. 9 and 11). 

3. The highest field intensities are seen in the resection cavity and residual tumor 
when the burrholes and electrode array are placed as close as possible above the 
mentioned tissue. This is illustrated when comparing mean values of the optimal 
electrode array placement of 60 degrees with a suboptimal placement of 150 
degrees (Figs. 10 and 13). 

4. The lowest electric field values in the residual tumor tissue and resection cavity 
were observed for positions 0 and 180 degrees (Figs. 11, 12, and 13). The mean 
range for those two electrode positions was for all tissues and SR-surgery posi- 
tions between 77.93 — 101.96 v/m. This indicates that electrode array positions 
that are far away or parallel with the burrholes should not be used since the direc- 
tion of the current does not pass as strongly through the region of interest, i.e., 
residual tumor and resection cavity. 

5. Placing the burrholes directly above the resection cavity increases the field 
strength in the residual tumor and resection cavity by 60—80% for most electrode 
array positions. The highest increase in the residual tumor was observed with the 
burrholes located at position (d) “up and back" and electrode placement 135-150 
degrees. This aligns with the burrholes and electrode array being close to directly 
on top over the residual tumor but not the resection cavity (Fig. 12). 


4 Summary and Discussion 


The five burrhole SR-surgery configuration was chosen for OptimalTTF-2 trial, 
where the central burrhole should be placed above the center of the residual tumor/ 
resection cavity. 

After the first trial patient underwent SR-surgery, the presented realistic compu- 
tational head model study was performed, mimicking the patient's tumor resection 
cavity, residual tumor, and location. The exact burrhole location and size were mim- 
icked from the postoperative CT scan of the head. 

The field strength was calculated with and without SR-surgery including elec- 
trode array rotation around the craniocaudal axis. The results indicated this 
SR-surgery configuration significantly increases the electric field strength in the 
residual tumor and resection cavity when the electrodes are placed near or directly 
above the burrholes. 

Finally, the SR-surgery configuration was moved to three alternative locations to 
examine if the central burrhole needed to be placed above the resection cavity as 
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—- Without burrholes 

—— With burrholes on top 

—— With burrholes moved back 
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Fig. 9 Cumulative distributions of electrical field strengths were obtained with burrholes directly 
above the resection cavity (red lines), the three alternative placements as shown in Fig. 8, and 
without burrholes (dotted lines). The cumulative distributions in the y-axis are given as the percent- 
ages of tissue (white matter, gray matter, resection cavity, and residual tumor) exposed to field 
strengths above the corresponding values in the x-axis. The position of the electrode array is 
selected to compare the optimal 60-degree placement and suboptimal 150-degree placement 
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Fig. 10 Axial view of 
TTFields intensity 
distributions in a head 
model that mimicked a 
trial patient’s postoperative 
MRI. The approximate 
positions of electrodes are 
shown as the surrounding 
patches. The left and right 
columns give the field 
distributions without and 
with burrholes, 
respectively. This show the 
field intensity distributions 
induced by electrode array 
positions: 0°, 45°, 90°, 
135°, and 180° (their 
corresponding electrode 
montages are given in 

Fig. 7). The field strength 
increases significantly 
when the electrodes are 
placed close to or directly 
above the burrholes 
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With burrholes on top 

With burrholes moved back 

With burrholes moved up 

With burrholes moved up and back 


Fig. 11 Statistical analysis of TTFields intensities with burrholes directly above the resection cav- 
ity (red lines), the three alternative positions as shown in Fig. 8 and without burrholes (dotted 
lines). The three columns represent the mean, median, and peak values of TTFields intensities for 
the four tissue types (WM, GM, residual tumor, and resection cavity). The x-axis shows the degree 
of rotation as seen in Fig. 7 and the y-axis shows the field strength in v/m. The peak values are 
defined as the 9946 percentile of the field intensities. The median and peak fields follow similar 
patterns for each tissue The TTFields intensities of WM and GM are similar with or without bur- 
rholes. The residual tumor and resection cavity field strength increase with burrholes, especially 
when they are close to or directly on top of the burrholes (corresponding to 60-150 degrees). 
Within this interval, the exact location of the electrodes can largely affect the peak values of the 
field intensities. The medians of field intensities are less sensitive to the electrodes' locations if 
they are set close to the burr holes 


initiated planned. For each alternative location, the rotation of the electrode array 
placement around the craniocaudal axis was performed. 

This concluded that the central burrhole and the entire SR-surgery configuration 
should be placed directly above the resection cavity and residual tumor since the 
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Without burrholes 

With burrholes on top 

With burrholes moved back 

With burrholes moved up 

With burrholes moved up and back 


Fig. 12 The mean and median E field intensity for WM, GM, resection cavity, and residual tumor. 
Without burrholes are shown as the dotted line and it can be seen at the level of the x-axis. X-axis 
also shows the degree of electrode array location. The y-axis represents the % increase in field 
intensity when SR-surgery is performed. The burrholes directly on top (red lines) and their alterna- 
tive locations 
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Fig. 13 Computational head model based on a healthy MRI with manually implemented features 
that mimicked a real patient’s postoperative MRI. Burrholes were manually added based on the 
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burrholes significantly increase the field strength of the pathological tissue directly 
underneath. 

These results that SR-surgery can improve the electric field and that it is possible 
to optimize the electrode array placement aligns well with previous research on this 
topic [5, 9, 10]. 

These results will form the basis for the standard operating procedure for 
SR-surgery configuration, its location, and electrode placement for each patient 
individually in the OptimalTTF-2 trial. 

There are a number of limitations to the work presented. This study was based on 
one computational head model study. The results and conclusion are therefore based 
on a limited amount of data that might not be valid if additional simulations were 
made. The residual tumor and resection cavity used in the head model might not 
accurately reflect the complex and heterogeneous morphology and location of all 
glioblastoma tumors. 

The five burrhole configuration was not tested in detail. Korshgj et al. conclude 
that several small burrholes are more effective than one big craniectomy, however, 
it remains unanswered what is the optimal amount of burrholes regarding field 
strength and without compromising patient safety. 

Furthermore, the results that the burrhole location should be placed above the 
pathological tissue seem intuitive. In scenarios where you only have a resection cav- 
ity or acomplete tumor intact, the chose where to place the burrholes will be simple. 
The burrhole placement choice is more uncertain when the choice is between a large 
residual tumor in the far periphery of the resection cavity. 

Therefore, future studies are needed to validate our findings in a larger number 
of heterogeneous models. We plan on examining the field strength for all the trial 
patients with and without SR-surgery, with realistic head modeling based on their 
MRI, CT, electrode placement, TTFields therapy data and compare with PFS, OS, 
and topographical areas of recurrence seen on the patients control MRI. 

Another point worth addressing is that the burrholes conductivity values were set 
to be the same as CSF, which may not accurately reflect reality, in which scar tissue 
may also constitute a part of this volume. This could give misleading results for 


< 


Fig. 13 (continued) postoperative CT-scan as shown in Fig. 4. The top three views (axial, sagittal) 
show a “funnel” or surgical access way and a resection cavity represented in blue. The tumor rem- 
nant is shown in yellow. The highest field strength was observed with electrode placement of 60 
degrees. This was compared with a suboptimal electrode placement of 150-degrees. The observa- 
tion is that the electric field strength is significantly enhanced in the tumor and resection cavity 
when electrode arrays are placed near or on the SR-surgery. Values: 60-degree electrode array 
placement. Burrholes directive above: Resection cavity 186.73 v/m, residual tumor 175.07 v/m, 
white matter 118.95 v/m, gray matter 91.96 v/m. Burrholes furthest away in the “back and up” 
position: Resection cavity 136.01 v/m, residual tumor 123.05 v/m, white matter 113.13 v/m, and 
gray matter 85.14 v/m. 150-degree electrode array placement. Burrholes directly above: 
Resection cavity 126.49 v/m, residual tumor 118.99 v/m, white matter 105.11 v/m, and gray matter 
83.25 v/m. Burrholes furthest away in the “back and up” position: Resection cavity 183.97 v/m, 
residual tumor 115.6 v/m, white matter 199.33 v/m, and gray matter 92.93 v/m 
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patients wearing TTFields therapy for months and years since their burrholes are 
covered in scar tissue. 
More research is warranted to address several of the above-mentioned issues. 


5 Conclusion 


In this chapter, we continue to explore the novel idea of skullremodeling surgery to 
optimize TTFields therapy and the importance of electrode array placement. The 
concept of SR-surgery (craniectomy, burrholes, and skull thinning) and optimal 
electrode array placement were examined in a phase 1 pilot trial, OptimalTTF-1 that 
concluded skullremodeling surgery to be safe and indicated an increased survival 
for first recurrence glioblastoma. 

The additional knowledge presented in this chapter is how skullremodeling sur- 
gery was attempted to be standardized to one SR-surgery configuration, it's place- 
ment and electrode array placement for the entire subsequent multi-center efficacy 
phase 2 trial “OptimalTTF-2”. 

Validating one SR-surgery configuration as the optimal would have a significant 
clinical impact since it would ensure uniformity of SR-surgery across the multi- 
center trial and thus increase the quality of the trial. 

A five burrhole skullremodeling surgery was proposed and with one realistic 
computational head modeling calculated to have a significant increase effect on the 
residual tumor and resection cavity region without increasing the field strength for 
grey- and white matter. The realistic head model was recreated using patient-specific 
MRI and CT data. There was a significant increase in electric field strength when 
the electrodes were placed near the skullremodeling surgery. The skullremodeling 
surgery technical aspect was described. The location of skullremodeling surgery 
mattered and the most optimal location was with the central burrhole placed above 
the residual tumor or resection cavity. The data generated from these simulations 
will be used for the standard operating procedure for skullremodeling surgery and 
TTFields electrode array placement, however, the results are based on limited data. 
Therefore, more research is planned by validating this approach from all individual 
patient data gathered from the trial as it progresses and depending on these results 
the standard operating procedure should be updated accordingly. 
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Part II 

Low Frequency Electromagnetic Modeling 
and Experiment: Neural Stimulation in 
Gradient Coils 


Peripheral Nerve Stimulation (PNS) A 
Analysis of MRI Head Gradient Coils gese 
with Human Body Models 


Yihe Hua, Desmond T. B. Yeo, and Thomas K. F. Foo 


1 Introduction 


With the recent advancements in high-performance head gradient coil technology 
for magnetic resonance imaging (MRI), the assessment of peripheral nerve stimula- 
tion (PNS) has become increasingly important. In this chapter, we first briefly 
describe the basic concepts of PNS and gradient coils in an MRI system. We then 
describe and compare the significance of a head-only MRI system with respect to a 
whole-body system. We will also provide a general outline of the electromagnetic 
simulation, the neurodynamic simulation and the gradient coil design processes. In 
the Methods section, we deliberate the shortcomings of a widely used human body 
model in the head gradient coil PNS analysis and describe the modifications made 
to two human body models to improve the correlation of simulation and measure- 
ment results from two head gradient coils. After that, we apply these new human 
body models to analyze three folded and non-folded gradient coils and reveal the 
relationship between the eddy current flow in the human body and the gradient coil 
wire pattern and its impact on PNS. We also discuss the connection between PNS 
and concomitant fields, and the effects that a human body model with simplified 
tissue properties has on the PNS calculations. Finally, we will discuss the shortcom- 
ings of this study and future work and provide our conclusions. 
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1.1 MRI Gradient Coil and E-Field 


An MRI scanner consists of several key hardware components, including the mag- 
net, the gradient coils and the RF coils. In a typical MRI scanner, the magnet 
provides a homogenous magnetic field, usually 1.5 Tesla or 3 Tesla, in a field of 
view (FOV) that spans 45 ~ 50 cm in diameter [1]. The gradient coil system has 
3 sub-coils: x, y and z. Each coil provides a linear B, field along one axis, such that 
B- = Gk (k =x, y, z) [2], where G, is the gradient strength of each sub coil. The RF 
transmit coil transmits RF waves at the Larmor frequency of the hydrogen nucleus 
to excite proton spins [2] while the RF receive coil measures the resultant MR signal 
from these spins. 

The resonance frequency o = yB, is proportional to the applied static magnetic 
field, where y is the gyro-magnetic ratio. The locations of the spins in space can be 
encoded by the B. field | using a combination of the main static field (Bọ) and the 


gradient by B, = B, + E G,k. In a modern MRI system, different gradient wave- 


forms are applied to each gradient sub-coil to spatially encode spins such that they 
have spatially varying frequencies and phases. The image of the subject can then be 
reconstructed from the received RF signals by utilizing image reconstruction algo- 
rithms [3]. 

The alternating magnetic field B generated by a gradient coil can induce an 
electric field E in our body, which can be calculated from the vector potential A 


E = - I ae 
and the scalar potential p as E = —0A/0t — VQ , with A(r) - a rr -F |d (F), 
T 


and t, J, dl, F, T being time, current in the gradient coil, differential element in 
length direction, source vector and observation vector, respectively. The quantity q 
in the body can be obtained from V -oVo+V -00A/60t 20 with j 2oE and the 
low-frequency approximation V - j 20, where c and j are electric conductivity and 
current density, respectively. Together with the boundary condition Lag +n a 0 

n t 
from fi-J =0,/ being the unit normal vector at the external surface, q can be finally 
computed with numerical methods such as Finite Element Method (FEM) for com- 
plicated human body geometry with heterogenous electrical conductivity for differ- 
ent tissues [4]. One of the main safety-related bioeffects induced by this E-field is 
the stimulation of our peripheral nerves, which may cause tingling or even painful 
muscle contractions in severe cases. Hence, PNS is an important safety topic in 
gradient coil design. All gradient EM simulations in this study were performed 
using Sim4life Ver.6.0 (Zurich Med Tech, Zurich, Switzerland). 
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1.2 Head vs. Whole-Body Gradient 


A very rough approximation of the E-field can give us some sense of the amplitude 
of the E-field. From Vx E -—0B/0t we get $ E-di = [Bas , where / is the 
t 


A 
integration path and the A is the area surrounded by the path, then [5] 
E(r)- — (1) 


It depends on the rate of change of the B-field and the volume (presented by r) to 
provide such a magnetic field. 

A gradient coil is composed of many turns of wire as shown in Fig. 6 of reference 
[6]. Therefore, the gradient coil has a non-negligible inductance, which indicates 
that the current through it cannot change instantly. This means that 0B/0t and the 
gradient slew rate (SR), denoted by dG/ot (subscript of G is removed later without 
loss of generality), cannot be infinite. A clinical whole-body gradient coil typically 
has a large FOV, resulting in a large-inductance coil that cannot support gradient 
waveforms with high slew rates. On the other hand, ultra-high performance head- 
only gradient coils designed for microstructure and functional imaging of the brain 
need a much smaller FOV that encompasses the head but require fast-changing 
gradient fields and high gradient waveform amplitudes. More importantly, head gra- 
dient coils with a smaller FOV will generate less E-field and thus have lower 
PNS risks. 

We recently developed two high-performance head gradient systems HG2 [7] 
and MAGNUS [8] for research purposes. The former is designed for a compact 3 T 
head MRI scanner and the latter is a head gradient that can be inserted into a clinical 
3 T whole-body scanner. Compared to typical clinical whole-body MRI scanners, 
whose gradient strength and slew rate are usually around 30 ~ 50 mT/m and 
100 ~ 200 T/m/s, HG2 and MAGNUS achieve 85 mT/m, 700 T/m/s and 200mT/m, 
500 T/m/s, respectively. In this study, these two coils are used in PNS model verifi- 
cation since we have measured data from volunteer scans. 

Figure 1 shows the E-field intensity of MAGNUS and a typical whole-body gra- 
dient coil with a head scan position in the Yoon-sun human body model (v4.0b03, 
IT'IS Foundation, Zurich, Switzerland). It can be clearly seen that the high | E| 
-field region of the whole-body gradient is much larger than that of the head gradi- 
ent. The high field regions by MAGNUS X and Y are mainly in the head. 


1.3 Peripheral Nerve Stimulation 


Nerve stimulation means an action potential is initiated and propagates through the 
nerve tissue. The theory is that the opening of the ion channels within the membrane 
of a neuron requires work from a mechanical impulse, which is further produced by 
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Fig. 1 E-field intensity of MAGNUS and common whole-body gradient coil with a head scan 
position in Yoon-sun human body model (SR is normalized to 100 T/m/s) 


an electric field through time. The fundamental law for nerve stimulation by an 
external electric field can be expressed as [5]. 


[E(r)atz E, (v. +7) (2) 


where E, is the rheobase, the minimum E-field to activate the nerve while t, is the 
chronaxie, the duration with which the threshold is exactly twice of the rheobase. If 
the external field is less than the rheobase E,, no stimulation can happen even if the 
field is applied for an infinite time. 

When the electric field is induced by the magnetic field, it is possible to use the 
magnetic field time derivative to describe the fundamental law. Since G and the SR 
dG/ot, are two important performance parameters for every gradient coil, the equa- 
tion can be further derived as a relationship between gradient strength, slew rate and 
the gradient pulse duration c [9]. 


AG(T,F)2 AG,, (Fr) G, (F)t (3) 


Because the E-field from the gradient coil for a certain scan position is diverse 
inside of the human body, the PNS response from each nerve is also different. The 
nerve that is easiest to be stimulated gives the threshold under this certain excitation 
condition. 
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1.4 Neurodynamic Simulation 


To study the action potential initiation and propagation in detail, a closer look at the 
behavior of the neuron membrane is needed. The sodium and potassium gates, etc., 
on the neuron membrane, can be simplified as a set of voltage-controlled conduc- 
tance, while other characteristics of the neuron can also be described as a set of 
distributed capacitance and conductance. A widely accepted computation model of 
myelinated mammalian nerve fibres is the MRG model [10], which utilizes double- 
cable RC-circuits to describe the Ranvier node, paranodal, internodal section of the 
axon and the myelin sheath of the axon. As a result, the action potential behavior 
can be simplified and described by the cable theory equation [11]. 


C, » + Gn Vi m G, [a = 2V, * Via) = G, (Vins E 2V,, T Va) (4) 


where Gn, Cm and G, are membrane conductance, membrane capacitance and axial 
internodal conductance, respectively. G,, is further controlled by membrane poten- 
tial reflecting the transient on-off status and the current-conducting ability of the 
sodium channel, potassium channel and other paths on the membrane. V,(t) and 
V.(t) are membrane potential response and applied potential by the external E-field, 
respectively (time dependency is omitted for simplicity in the equation). The excita- 
tion item (driving function) on the right side of the equation is about the second 
order difference of the electric potential, but actually relates the first order differ- 
ence of the tangential E field along the nerve, where n is the index of the Ranvier 
node. With the E-field calculated from the electromagnetic simulation in the previ- 
ous step, the voltage drop V, along the nerve trajectory can be obtained by 


V, = fE -dl . After multiplying with the time domain modulation function, a(t) 


B-field by Gradient Coil 
E-field a(t) 


L 


Fig. 2 Curves to represent the relationship between the B-field waveform and E-field waveform 
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Table 1 Parameters used in neuron simulation 


Neuron model | MRG 
Diameter(um) 16 (motor) 
12 (sensory) 
B ramp time c (ms) | 0.2.0.5 
B plateau time (ms) | 1 
Train length of pulses | 16 


(E-field waveform through time, refer to Fig. 2), and an additional scaling coef- 
ficient T, 


VO =V, a(t): T (5) 


is finally fed into the Eq. (4) to calculate the membrane potential response for each 
nerve trajectory. The solver will iteratively change the scaling factor T to find the 
smallest field intensity to initiate an action potential. This is called the titration pro- 
cess. In this study, we use the NEURON module (Yale University, New Haven, CT) 
that is integrated in Sim4life for PNS calculations. Thus, AG can be calculated for 
every c. Since the relationship between AG and 7 is linear, only two values of t are 
used in the simulation. Table 1 shows the parameters used for neuron simulation in 
Sim4Life. 


1.5 Gradient Coil Design 


Modern gradient coil design can be regarded as an optimization problem that seeks 
the optimal current distribution on a pre-defined design surface that produces a tar- 
get field distribution while subjected to a set of optimization constraints [12, 13]. 
The current density can be decomposed into basis functions with unknown 
coefficients. 


J 2Xcf.i-l..n (6) 


and the interested physical quantities such as magnetic energy and the B, field in 
FOV can be expressed as quadratic form or linear form of the unknown coefficients, 
respectively 


1 1 J(F)-J(F ) „1 
reir! xg ard’ = c Le (7) 


g (y) - He peer] "T m 
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for i, j = 1, ..., n. Finally, the problem usually could be solved by a quadratic pro- 
gramming algorithm: 
To find: c = [c), C2, ..., Cal” 


NU 1 : ; 
Minimizing: ae c + other quadratic forms of the unknown coefficients... 


Gradient constraint: G.x-e <B.c<G,x+e for each point in X coil FOV, e is 
the tolerance. 

Other constraints such as force, torque, eddy current and fringe field having lin- 
ear form of the unknown coefficients can also be included. 


2 Methods 


2.1 Augmentation of Nerves in Yoon-Sun and Duke Models 


We used Yoon-sun female model v4.0b03 and Duke male model v3.1b01 from IT'IS 
foundation in this study. The Yoon-sun model is reconstructed from high-resolution 
cryosection images and it has a relatively complete nerve trajectory atlas for further 
neurodynamic simulation while the Duke model is based on MRI data of a volunteer 
without any nerve trajectories. In our previous study of PNS response for MAGNUS 
[14], we found that the calculated PNS thresholds for the X and Y coils are higher 
than those in the measured data. We hypothesized that this was because the existing 
nerve trajectories in the head of the Yoon-sun model are located mainly in the deeper 
recesses of the skull, which are in the low gradient field regions. The Yoon-sun 
model lacks extracranial nerves that are located in the high gradient field region. 
To address this, we added extracranial nerves and superficial nerves in the neck 
and shoulder region to both Yoon-sun and Duke based on our knowledge of the high 
field area of the head gradient coils and the anatomy of the body (Fig. 3). The reason 
we only added a portion of the nerves is that the human body and the X-coil are 
mainly left-right symmetric, and it is time-consuming to add anatomically realistic 
nerve trajectories into the existing human body model. Table 2 summarizes the 
nerve branches that have been added to the human body model. The nerve branches 
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Original 929 traj. 


Fig. 3 (a) Yoon-sun: Original nerve trajectories (White), 86 nerve trajectories (Yellow) added and 
(b) Duke: 83 nerve trajectories added 


Table 2 Added nerve branch names 


Nerve name 
Buccal branch of left facial 


Temporal branch of left facial 


Lateral branch of left supraorbital 
Medial branch of left supraorbital 
Left supratrochlear 


Left infratrochlear 


Left recurrent laryngeal 


Right recurrent laryngeal 


Right dorsal scapular 


Right medial supraclavicular 


Right accessory 


Cervical branch of right facial 


selected in Duke differ slightly from Yoon-sun based on our early calculations on 
the latter by omitting nerves that have a low likelihood of being stimulated. 


2.2 Non-folded and Folded Gradient Coil Design 


We used the gradient specifications (Table 3) from Tang's paper [12] to design three 
different gradient coils for this study (Fig. 4). Coil A is a non-folded design, in 
which the primary layer is separated from the shield layer. The second design is 
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Table 3 Parameters for designing folded and non-folded gradient coils 


Lprimary (cm) 50.4 Dyoy (cm) 24 

Lsnieta (cm) 72.6 Offset of FOV center to coil geometry center (cm) —9.2 
Lsm) ———— (137. (Curn(A) —— A 

D primary (cm) 33.6 Gradient strength (mT/m) 60 

Dyhieia (cm) 55.8 € (max field error in FOV) 0.1 

Drs (cm) 90.6 A (max field error by eddy current) 0.002 


Fig. 4 Three X coils designed for this study. (a) A Non-folded design; (b) a folded design; and (c) 
a folded design with turn number limit in primary/shield connecting region 


folded at the patient end, allowing more turns close to FOV which makes it more 
efficient. However, the wire turns that are close to the FOV are also close to the 
human body, which increases the E-field exposure. The third coil is also folded, but 
the coil turn number in the primary/shield connecting region is reduced to about half 
of the second coil. The idea is to seek a balance between the first and second design, 
and to examine how the wire pattern strategy affects the PNS in the human body. 

It needs to be noted that PNS constraints were not applied in the gradient coil 
optimization process in this study. Relatively mature algorithms have been devel- 
oped during recent years and an interested audience could read Davids’ work [15] 
for further information. In this work, however, we mainly focus on the mechanism 
of how the wire pattern changes can affect the E-field, and thus, PNS. 


3 Results and Discussions 


3.1 Model Verification on MAGNUS and HG2 X Coils 


The simulated PNS threshold by MAGNUS and HG2 X coils are within the observed 
range and have shown much better agreement with the measurement than before 
(Fig. 5a) [14]. Figure 6 shows the most sensitive places for PNS in Yoon-sun and 
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PNS Model Validation PNS results in Folded/Non-Folded Coils 
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Fig. 5 (a) PNS simulation results compared to measurements on MAGNUS and HG2; (b) PNS 
simulation results in three folded and non-folded gradient coils 


Fig. 6 PNS initiation places on nerve trajectories stimulated by MAGNUS (a, c) and HG2 (b, d) 
are represented by small balls. (a, b) are with Yoon-sun models (c, d) are with Duke model. The 
bigger and the redder of the ball means the location is more sensitive 


Duke. Overall, MAGNUS and HG2 X coils have very similar PNS performances. 
PNS was observed to be strongest at the nasion/glabella region and also occurs with 
high probabilities in regions close to the eye, forehead and cheeks [16, 17]. These 
places are obviously different from those in [14] since there were no such nerve 
trajectories in the original Yoon-sun models. PNS threshold is lower in Duke than in 
Yoon-sun, mainly because the size of Duke is larger than Yoon-sun, and as such, the 
nerves in Duke can extend to higher E-field regions than in Yoon-sun. All these 
results are consistent with observations from volunteer scans that were reported in 
our previous work [16]. 

We also observed that simulated thresholds are lower than the measurements. It 
is probably because, in simulations, the action potentials are initiated mostly at the 
ends of the nerve trajectories, which seem to be induced by the rather abrupt termi- 
nation of the nerve trajectories that is unreal. In reality, these nerves may extend 
longer and become narrower in diameter, having telodendria or dendrites at the 
ends. We have reported this effect in our previous study [14], but for comparison 
purposes, the simulation result is still close to measurement. 
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3.2 B,and |E | Field in Free Space and Human Body 
for Folded/Non-folded X Coils 


The Bz and E-field distributions in free space in the three gradient coil designs are 
plotted in Fig. 7. While the B-field will be largely unperturbed by the presence of a 
human body model, the E-field distribution could be very different when a body 
model is introduced into the simulations. It could be seen that for the three designs, 
although the Bz fields are similar in FOV, the E-field in free space and the human 
body could be different. Coil A generates a relatively large E-field in the head 
region, while coil B results in a larger E-field in the neck and shoulders. Coil C low- 
ers the E-field in the body from coil B by limiting the coil turns in the fold region. 
We should also note here that due to interference of the shoulder of Yoon-sun and 
the gradient coil, Yoon-sun has to be shifted 4 cm inferior to Duke in these three 
coils. This is another reason why augmentation of nerves was introduced into both 
Yoon-sun and Duke models. 


Z [m] 


Z [m] 


X [m] X [m] X [m] 


Fig.7 Field of the folded and non-folded gradient coils in free space. Upper row: Bz in XZ plane 
(scaled to 10mT/m); lower row: | E| in XZ plane (scaled to 100 T/m/s). Red arrows show the dif- 
ferent field distributions in three coils 
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3.3 PNS Simulation Results for Folded/Non-folded X Coils 


While Figure 8 provides the E-field intensity, Figure 9 shows the most PNS-sensitive 
locations in Yoon-sun and Duke in the three newly designed gradient coils. For Coil 
A, the action potentials are initiated mostly in the forehead region in Yoon-sun. For 
the same coil, the action potentials are formed in both the forehead and glabellar 
regions in Duke. For coil B, the strongest PNS is observed in the neck or shoulder 
region. For Coil C, the most PNS-sensitive regions have a situation between coil A 
and coil B. Coil C yields the highest threshold so the lowest risk in the Yoon-sun 
model but has similar performance to coil B for Duke. For coil C, although the 
E-field in the body is higher than the E-field in the head (Fig. 8), the PNS threshold 
in the body is not necessarily lower. This is because, from the nerve cable equation, 
it 1s the first order difference of E-field intensity along the nerve that determines 
PNS, not the E-field intensity itself. But overall, PNS-sensitive locations are corre- 
lated to the high field region of the gradient coil. The difference in PNS results 
between the models can be explained by the difference in the size of the subjects 
and their relative positions in the coil. Figure 5b summarizes the PNS thresholds of 
the three coils. 


Coil A Coil B Coil C 


21 0 01 02 22 21 0 01 02 22 o1 ] 01 02 


Fig. 8 |E| displayed on XZ plane, upper row: Yoon-sun, glabella@-60 mm from gradient coil 
ISO (due to the coil/shoulder interference); lower row: glabella@-20 mm from gradient coil ISO 
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CoilA Coil B Coil C 


Yoon-sun 


Duke 


Fig.9 Locations of PNS initiation in nerve trajectories stimulated by three X coils are represented 
by small balls. White arrows mark the most sensitive regions in case difficult to read 


3.4 E-Field Streamlines 


We recommend using E-field streamline plot to further reveal the relationship 
between the gradient coil wire pattern and PNS as what is shown in Fig. 10. E-field 
follows the eddy current J in the human body since J =o £, so that using E-field 
streamline is a better way to represent how the eddy current flows in the human 
body than E-field intensity map. It can be seen that eddy current should follow the 
wire pattern of the gradient coil and is also affected by the geometry and tissue 
properties of the human body. 

For Coil A, since most of the conductor turns are at the upper head region, higher 
intensity E-fields and thus eddy currents, will flow through this region. On the other 
hand, since Coil B has many wire turns in the lower part of the coil, the eddy current 
loops concentrate on the neck and shoulders. This also leads to reduced eddy cur- 
rents in the forehead region. Coil C yields a situation in between Coils A and B. It 
is clear that the E-fields that form the eddy current determine the PNS. Thus, to 
design the gradient coil with PNS constraints, one would seek to shape the eddy 
current paths in the human body to avoid producing high E-fields (driving function) 
in regions where peripheral nerves are located. While previous studies propose 
algorithms to constrain PNS in the gradient coil design process, here, we provide a 
more intuitive explanation of how they work. 
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Fig. 10 E field streamlines in Yoon-sun model by three new designed coils: (a) coil A, (b) coil 
B, (c) coil C rendered with overlapping wire pattern of the coils (SR scaled to 100 T/m/s) 
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Fig. 11 Concomitant Bx field rendered together with Bz field for three new designed coils 


3.5 Concomitant Field 


Conventionally, only B, is controlled in gradient coil design since only B, is involved 
in image encoding. Previous studies [18, 19] showed it is possible to reduce the PNS 
risk by manipulating the concomitant field of the gradient coil. The concomitant 
field exists because the magnetic field is harmonic at places of no source, and it is 
impossible to obtain a linearly changing B, without involving changing B, and/or B, 
[20]. Figure 11 shows the concomitant B, fields together with the B, fields for the 
three X coils. We can see that while B. fields of these coils are identical in FOV, the 
B, fields are different. The correlation between | E| and | B | can again be roughly 
explained by Eq. (1) and be verified with Figs. 11 and 12. | E| field can be further 
manipulated by involving an additional uniform field coil [19] since essentially G is 
used in MR imaging process. Basically, the concomitant field can be changed by the 
changing of the wire pattern, as the same as that the eddy current loop in the human 
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Fig. 12 Free space |B field of three folded and non-folded coils in XZ plane (SR scaled to 
10 mT/m). Coil A provides relatively higher | B| values in the upper FOV region and Coil B and 
Coil C provide relatively higher | B| values in the lower part (neck and shoulder region) 


body can be changed by the wire pattern. But it should be noted that directly con- 
trolling the concomitant field does not involve the information of the human body, 
which includes the geometry and the tissue properties. 


3.6 Effects of Homogeneous/Simplified Tissue Properties 


The homogeneous human body model is attractive because it may simplify numeri- 
cal computations and enable the utilization of the boundary element method (BEM), 
which might be faster than the FEM. But BEM only has a computational advantage 
when the boundaries or numbers of tissue types are relatively small. In the IEC 
60601-2-33 standard, it has been recommended to use a homogeneous cylinder for 
a rough estimation of the | E | field. Although 0.2 S/m is mentioned in this standard, 
which is roughly the average electrical conductivity in the human body, the E-field 
will not be affected by any electric conductivity value for the homogeneous model. 
Another situation is to use the human body model with as few tissue properties as 
possible. This is interested because segmentations for different tissues are very time 
expensive. The creation and use of subject-specific human body models for accurate 
PNS estimation in high-performance gradient systems are gaining more attention, 
which parallels the need for similar models for local specific absorption rate (SAR) 
estimation in ultra-high field MRI. The original Yoon-sun and Duke model are com- 
posed of over 70 tissues of different electrical conductivities. We conducted several 
simulation cases to study the impact on PNS results by simplifying the tissues, 
including (a) a homogenous human body model; (b) a model with only 3 tissue 
properties (fat, muscle and bone, all the other tissues are treated as muscle); (c) a 
model with 6 tissue properties (lung, liver and skin are further distinguished accord- 
ing to case (b)) and (d) a model with 8 tissue properties (averaged electric 
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conductivities of tissues inside of the skull and tissues inside of the belly are further 
added in simulation to represent tissue inside of skull and tissue inside of belly, 
respectively). 

The PNS responses of simplified tissue models are compared with the original 
heterogeneous Yoon-sun model on MAGNUS X and Z coils. It is noted that it is 
beneficial to retain accurate tissue model properties for tissues in closest proximity 
to nerve tissue. Tissues that are far away from the nerves can be assigned tissue 
properties that are averaged across the local tissue types. This is a compromise that 
reduces computation time while increasing the likelihood that the E field distribu- 
tion in closest proximity to nerve tissues conform closely to that in a fully seg- 
mented human body model. Table 4 lists the conductivities used in the simplified 
models. To compare the PNS response from different cases, we first convert the 
PNS response lines of different nerve trajectories in t — AG plane to be points in 
AG nin — SR plane. Then we examine how the responses of the nerves shifted in the 
AG, — SR plane due to the property simplification. We can calculate the pseudo- 
distance between the corresponding points to quantify the difference but only the 
most sensitive nerves (for example 20) matter. We can also count how many trajec- 
tories remain in the list of the original 20 most sensitive nerve trajectories. As an 
example (Fig. 13 and Table 5), for the 8-tissue model, stimulated by MAGNUS X 
coil, the averaged distance is 13.61 for the 20 most sensitive nerves and 20 of them 
are remained as the most sensitive nerves in Yoon-sun, while these numbers for the 
homogenous model become 48.62 and 16, respectively. Table 5 summarizes all the 
cases performed for the Yoon-sun model. It shows that compared to the homoge- 
nous model, the 3-tissue model provides a closer result to the original heteroge- 
neous model and that the 6-tissue model is even closer to the latter. The 8-tissue 
model is a little better than the 6-tissue one because the averaged electric conductiv- 
ity for tissues in the skull and in the belly is close to that of the muscle. Another 
Observation is that fewer tissue properties are needed for the MAGNUS X coil than 
for the Z coil. This is because there are inherently fewer types of tissue outside of 
the skull than in the body. Overall, the result shows that the homogenous model 
could be used to give a rough estimation on the E-field but might be not accurate 
enough to estimate the PNS. Six tissue properties or more might provide a better 
estimation. 


Table4 Tissue conductivities used in simplified tissue models 


Tissue name o (S/m) 
Fat 0.057 
Muscle 0.355 
Bone 0.004 
Skin 0.170 
Lung 0.105 
Liver | 0.220 
Averaged tissue in the skull 0.312 
Averaged tissue in the belly | 0.310 
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Fig. 13 PNS pseudo-distance of the 20 most sensitive nerve trajectories between simplified tissue 
model with 8 tissues and the original heterogonous Yoon-sun model. (The number marked on each 
point is the ID assigned to the nerve trajectory. Different colors are for different nerve branches) 


Table 5 The averaged pseudo-distances and top sensitive nerve numbers (of the 20 most sensitive 
trajectories in the heterogonous model) of simplified tissue models 


Simplified tissue Averaged Top sensitive nerve no. 
model pseudo-distance (of 20) 

MAGNUS |8 tissues 13.61 20 

X 6 tissues 16.61 19 
3 tissues 32.68 17 
Homogeneous 48.62 16 

MAGNUS Z | 8 tissues 5.76 19 
6 tissues 5.9 19 
3 tissues 57.16 8 
Homogeneous 110.42 8 


4 Shortcomings and Future Works 


There are two shortcomings in this study. The first is that only a portion of the nerves 
has been reconstructed, limiting PNS analysis for X Coils. More nerves are needed to 
enable full X/Y/Z coil design and analysis. Moreover, the nerve trajectories’ geome- 
try might not be perfectly accurate anatomically since a detailed image-derived nerve 
model was not available. However, compared to the lack of nerve trajectories in the 
interested places in the original body models, these added trajectories provide richer 
information to capture PNS responses for gradient coil analysis and design. A similar 
situation has been reported in [17], where virtual human bodies models are used. The 
second shortcoming is that the PNS constraint is not yet integrated into the gradient 
coil design in this study. It is preferred to include this constraint in the future. 
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5 Conclusions 


In this work, extracranial nerves and upper body superficial nerves have been added 
into both Yoon-sun and Duke models. PNS analysis for gradient coil in the human 
body models can give results of reasonable accuracy with a relative complete nerve 
trajectory atlas in the places of interest. Three folded and non-folded head gradient 
coils of different wire patterns have been designed to reveal the relationship between 
the E-field/eddy current flow in human body and PNS. E-field streamline plots 
inside of the human body can provide more information than intensity projection 
plots. To manipulate PNS through concomitant field generally is to tune the E-field 
distribution through changing the B-field distribution but does not include the infor- 
mation of the human body. A homogenous model is not accurate enough for PNS 
threshold estimation, instead, six or more tissues are needed for simplified tis- 
sue models. 
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1 Introduction 


Transcranial Magnetic Stimulation (TMS) is a non-invasive and safe method for 
activation of cortical regions by delivering high amplitude, short current pulses into 
a coil adjacent to the subject’s scalp. This creates a strong time-varying magnetic 
field (~1—2 T), which in turn induces an E-field at the surface of the brain [1, 2]. The 
E-fields can depolarize/hyperpolarize the cell membrane to activate/inhibit neuronal 
populations [3]. TMS is approved by the FDA for treating neuropsychiatric disor- 
ders such as major depressive disorder (MDD) [4] and obsessive-compulsive disor- 
der (OCD) [5], with more clinical applications under investigation. While accurate 
targeting of a brain region with a single TMS coil is a highly useful and accurate 
method for various experimental paradigms, characterising the functional connec- 
tivity of adjacent areas of a cortical network is only achieved by dual or multichan- 
nel TMS systems [6, 7]. Additionally, multichannel TMS arrays allow electronic 
steering of the induced Electric field (E-field) without the need to physically move 
the coils [2, 8] resulting in the capability of rapidly shifting the target/focus. This 
can be accomplished by computationally determining the E-field distribution of 
each coil in the array individually and then combining them all together to 
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synthesize a desired E-field pattern to stimulate any given cortical target according 
to a user-specified cost-function [2]. 

In general, a major challenge in TMS targeting is how to position the coil with 
respect to subject's head for optimal stimulation of a desired brain region [9-11]. To 
address this, a neuronavigation system can be used in which an optical tracking 
device localizes the TMS coil with respect to the subject's head [12], ensuring con- 
sistent coil placement across multiple TMS sessions [13]. However, using a basic 
navigation system without anatomically realistic E-field modeling, the actual stimu- 
lation intensity at the intracranial target and the surrounding regions remains 
unknown. Therefore, it is critical to computationally estimate the distribution of the 
TMS-induced E-field intensity in order to delineate the affected brain regions in 
real-time [12]. In general, we need to consider several factors for precise computa- 
tional TMS navigation. First, the anatomy of subject's head must be accurately reg- 
istered to their MRI data, ensuring that systematic coil positioning errors are 
minimized [11]. Second, we need to rapidly and accurately calculate the E-field 
pattern of the TMS coil or coil array to create the best possible “match” with the 
target area [12]. The real-time E-field based computational neuronavigation is uti- 
lized in our multichannel TMS system within a slightly different way such that 
when the coil array position is fixed with respect to the subject's head, the E-fields 
need to be computed only once and different current amplitudes will be applied to 
each coil to synthesize a desired E-field pattern at the target [8]. However, if the 
subject's head moves between pulses, the currents need to be adjusted to compen- 
sate for this, requiring a "near real-time recalculation" of the induced E-fields for all 
the coils in the array. 

For interactive E-field based navigation of the TMS coil during a stimulation ses- 
sion, the E-field distributions need to be ideally computed and displayed within 
100 ms (frame rate of 10 Hz). While simple (spherical) models used in commercial 
neuronavigation systems allow fast E-field computations, they may give inaccurate 
results due to oversimplification of the volume conductor (head) model [14, 15]. On 
the other hand, while high resolution individualized head models used in Boundary 
Element Method (BEM) and Finite Element Method (FEM) reach higher spatial 
precision, they are computationally too sluggish for real-time navigation requiring 
- 10 seconds per solution [16-18]. We have previously addressed this speed vs. pre- 
cision dilemma by introducing a new method based on the concept of an individual- 
ized Magnetic Stimulation Profile (MSP) [19]. The real-time E-field estimation by 
MSP approach leverages pre-calculation of the E-fields from a set of dipoles placed 
around the head model. The pre-calculations are done with the Boundary Element 
Method accelerated by the Fast Multipole Method (BEM-FMM) [20, 21]. Since the 
total E-field of an arbitrary TMS coil only depends on the incident E-field and tissue 
conductivity boundaries [20], we can find the matching coefficients between the 
incident E-field of the coil and the dipole basis set approximation, and the total 
E-field of the coil will be obtained by applying these matching coefficients to the 
total E-fields of the dipoles [19]. 
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Here we used two 3-axis coils [8] to illustrate how different combinations of the 
array elements allow steering the E-field ‘hot spot’. Furthermore, our computational 
multichannel TMS neuronavigation system allows rendering the E-fields of indi- 
vidual 3-axis coils and synthesizing the ‘hot spot’ using the array approach with 
speed and accuracy suitable for human studies. To verify that the computational 
neuronavigation system is working as expected, we used two z-elements of our 2x3- 
axis coil array connected to individual stimulators to mimic a figure of eight coil to 
activate the motor cortex and compared the results with a commercially available 
TMS coil. 

This article is organized as follows: In Sect. 2 we describe the design of the 
3-axis coils and investigate the efficiency of these coils in a simple 2x3-axis array 
configuration by calculating the induced E-field using a spherical model. Section 3 
provides a brief description of TMS neuronavigation methods as well as their pres- 
ent limitations in practical applications. In Sect. 4 we briefly describe the methodol- 
ogy behind the MSP approach for real-time calculation of the E-fields and its 
integration with a commercial TMS navigation system. Finally, in Sect. 5, we dem- 
onstrate the efficiency of two z-elements in a 2x3-axis array in conjunction with the 
interactive navigation system to elicit Electromyography (EMG) responses in a 
healthy volunteer. 


2 Multichannel TMS Array Design Concept 


2.1 The 3-Axis Coils as the Building Blocks of a TMS Array 


The proposed modular 3-axis multichannel array allows efficient and safe stimula- 
tion of any cortical area with high degrees of freedom in shaping and adjusting the 
orientation of the E-field. Delivering specific current amplitudes for each coil 
enables simultaneous or sequential stimulation of multiple areas which can be uti- 
lized for investigation of causal relationships between cortical areas involved in 
various information processing tasks [22, 23]. 

Figure 1a shows the basic conceptual design of our envisioned 48-channel (16 x 
3) TMS array. In this array configuration that has been motivated by our overarching 
goal of an integrated TMS-compatible MRI acquisition system [24], there are 16 
locations for placement of 3-axis coils allowing whole-head coverage. The 3-axis 
coil consists of three orthogonal circular elements called X, Y and Z as shown in 
Fig. 1b, c. The goal of this coil design is to achieve high efficiency and focality 
while offering accurate spatial control of the E-field hot spots using simple coil 
modules. The purpose of adding the X and Y elements is to cover for the zero E-field 
at the center of the circular Z-element [8]. Additionally, the X and Y elements will 
provide more degrees of freedom to the overall field shaping capabilities of the 
multichannel array. Furthermore, the geometrical design of the X, Y and Z elements 
eliminates the coupling between them within each unit [8]. 
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A) 


C) X-element 


Y-element 


Z-element 


Fig. 1 (a) Multi-channel TMS coil array. (b) prototype of the first of its kind 3-axis coil. (c) Each 
3-axis coil consists of three orthogonal elements allowing separate and combined stimulations 


2.2 The Numerical Method for E-Field Computation 
of TMS Coils 


The calculation of the E-field induced by each coil element was done using the 
BEM approach accelerated by the Fast Multipole Method (BEM-FMM) [20]. Based 
on the Maxwell-Faraday law, the magnetic field of the TMS coil induces an elec- 
tric field: 


vena (1) 
ar 


where B is the magnetic field of the coil. The total electric field E can be written as: 
E=E™ pi - P vg Q) 
t 


The E"* = — dA/ot corresponds to the primary or ‘incident’ field induced by the 
current in the coil, and the E* = — V q is called the secondary field that arises from 
the differences in the electrical conductivities of various tissues inside the head [25]. 
The BEM-FMM operates with tissue conductivity boundaries where E"* causes 
accumulation of surface charges giving rise to the secondary field — V @. The 
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boundary integral equations stemming from the quasi-static current conservation 
conditions V - J = 0 are numerically solved by the Generalized Minimal Residual 
Method (GMRES) [26]. Once the solution converges and the charge distribution on 
the conductivity boundaries is known, the total E-field in any region of the 3D space 
can be obtained. 

The individual E-field patterns generated by each element were previously char- 
acterized [8] and here we explore different E-field patterns obtained by various 
combinations of the elements. The coil models were generated in MATLAB [21] 
(MathWorks, Inc., Natick, MA, USA), placed over a 1-layer spherical model and 
the E-fields were calculated on the inner sphere at approximately 2 cm distance 
from the center of the coils. Figure 2 shows the model configuration as well as 
several calculated E-field patterns based on the activated elements and current 


dt 

of Maximum Stimulator Output (MSO). However, for actual stimulation experiments, 
the current intensity delivered to each element can be optimized using Minimum 
Norm Estimate (MNE) method [27] in order to generate a desired E-field pattern. 

Figure 2 also shows the E-field pattern of a commercially available figure-of- 
eight coil (C-B60, MagVenture, Farum, Denmark) as a reference. A similar but 
somewhat more focal pattern can be generated by combining two of the z-elements 
(Fig. 2, b8). However, the C-B60 generates a fixed E-field pattern with maximum 
intensity at the intersection of the two wings and to stimulate a different region the 
coil must be physically moved. Generally, even a simple multichannel array 


polarities. The coil current rate of change (2) values were corresponding to 50% 


Fig. 2 (a) Two 3-axis coils positioned on a spherical head model. (b1-9) Corresponding E-field 
distributions of the two 3-axis coils with the activated elements shown in green. (c) A commer- 
cially available standard figure-of-eight TMS coil model (Mag Venture C-B60) and the correspond- 
ing E-field pattern (d) 
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comprising of two 3-axis coils provides high flexibility in E-field shaping while still 
offering stimulation intensities comparable to commercially available figure-of- 
eight coils. 


3 TMS Navigation Systems 


The rationale of using an image-guided TMS navigation system is to reduce the 
variability of TMS coil placement with respect to the subject's head and hence to 
improve the repeatability and accuracy of TMS targeting and dosing [28, 29]. TMS 
navigator systems typically use optical tracking devices to localize the coil with 
respect to the subject's head. Anatomical landmarks on the subject's head are identi- 
fied with corresponding points in the Magnetic Resonance Imaging (MRI) data of 
the subject, allowing co-registration of the TMS coil position/orientation with the 
individual anatomy. 

Here, we used a commercial TMS neuronavigation system (LOCALITE GmbH, 
Germany) that employs an optical tracking camera (Polaris Spectra, Northern 
Digital, Inc. Waterloo, Ontario) to capture the position of the reference and coil 
trackers and displays the information graphically within the neuronavigation soft- 
ware. Additionally, a specific stimulation target can be defined using the software. 
The interactive navigation tool allows for accurate positioning of the coil at the 
target during the stimulation experiment. The positioning of the coil with the neuro- 
navigation system is naturally prone to some errors such as the method used for 
optical tracking of the coils and the head as well as the registration accuracy of the 
subject's head to the MR data [30]. Moreover, quantitative high-resolution target- 
ing/dosing of TMS requires knowing the E-field intensity at the desired cortical 
location and surrounding regions [31] which is necessary especially for multichan- 
nel applications but not available on commercial navigation systems. Here, we lev- 
eraged our recently developed MSP approach [19] for near real-time TMS E-field 
calculation combined with a commercial navigation system to track the coil move- 
ments and render the induced E-field with a frame rate of 6 Hz. 


4 Fast E-Field Calculation with Dipole Based Magnetic 
Stimulation Profile Approach 


A detailed description of the MSP approach can be found in [19]. The computa- 
tional pipeline of the MSP approach consists of a pre-calculation and a real-time 
step. In the pre-calculation step, hundreds of magnetic dipoles are distributed uni- 
formly on the scalp model as shown in Fig. 3a, and the incident (field from the 
dipoles only) and total E-fields (fields from dipoles coil and charge accumulation at 
tissue conductivity boundaries) are computed using the BEM-FMM method [20]. 
This fundamental basis solution only needs to be calculated once per subject and 
can be subsequently used to quickly estimate the E-field created by any TMS coil. 
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Fig. 3 (a) Spatial configuration of an example basis set with 3 orthogonal magnetic dipoles posi- 
tioned at several hundreds of locations around the subject’s scalp surface for pre-calculation of the 
incident and total E-fields. (b) A model of a commercially available TMS coil (C-B60, MagVenture, 
Farum, Denmark) is placed at an arbitrary location and the dipole amplitudes (c) are optimized to 
match the incident E-field of the coil. (d, e) the incident and total E-field of the C-B60 coil calcu- 
lated by the BEM-FMM as the ground truth. (f) The incident E-field obtained from the dipoles 
basis set. (g) The matching coefficients obtained in (c), are applied to the dipole basis set total 
E-fields to approximate the total E-field of the coil 


The real-time step is based on fundamental physical principles of TMS-induced 
E-fields: 


* The total E-field of a TMS coil only depends on the incident field and the tissue 
conductivity boundary surfaces and the associated conductivities [1]. 

* The incident E-field created by the coil in free space can be re-computed very 
quickly by applying translations/rotations to the set of points where the Incident 
E-field was initially calculated on. 

* By projecting the incident field of the TMS coil onto the set of dipoles and 
obtaining a set of weights that provide the optimal *match', the same weights can 
be used to approximate the total E-field based on the principle of superposition 
[32,33]. 


Based on these principles, the pre-calculated E-fields of the dipole basis set can be 
utilized for a very fast estimation of the total E-field of an arbitrary TMS coil as 
shown in Fig. 3. This method allows for real-time (100 ms) computational perfor- 
mance and display of the TMS-induced E-field for interactive visualization of the 
stimulated cortical regions. 

Figure 4 shows an example of the real-time setup consisting of a MSP-based 
computational E-field engine integrated with a commercial navigator system 
(LOCALITE, Bonn, Germany), using a head phantom. The reference head tracker 
placed on the phantom allows for co-registration of the head shape with a 'synthe- 
sized MRI dataset’. The position of the coil with respect to phantom is recorded and 
displayed by the navigator software using the tracker attached to the coil and an 
optical position measurement camera. The position of the coil(s) (either the two 
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Fig. 4 The real-time E-field computational modeling setup consisting of a commercial neuronavi- 
gation system (LOCALITE, Bonn, Germany) interfaced with an external PC for MSP-based cal- 
culation and display of (a) a commercial figure-of-eight coil (C-B60) and (b) two custom-made 
3-axis coils (Tristan Technologies, San Diego, California) 


z-elements of the two 3-axis coils or the C-B60 coil) are then streamed into a second 
PC using the TCP/IP and JSON protocols implemented in the LOCALITE naviga- 
tor. Finally, with the coil position/orientation streamed continuously to MATLAB, 
we can calculate the total E-field and display the results on a cortical surface as 
shown in Fig. 4. In the current setup, we were able to freely move the coil around 
the phantom and calculate/visualize the total E-field within a frame rate of 3 Hz and 
6 Hz for the two z-elements of the two 3-axis coils and the C-B60 coil, respectively, 
using MATLAB 2021a on an Intel Xeon(R) Gold 6226R PC with 192 GB of memory. 


5 Motor Cortex Stimulation Experiment with Z-Elements 
of Two 3-Axis Coils 


We recruited one healthy male volunteer to evaluate the performance of the simul- 
taneous activation of currents of opposite polarity in the Z-elements of two 3-axis 
coils in a motor cortex stimulation experiment. Informed consent was obtained from 
the participant in accordance with the study protocol that was approved by the local 
IRB at Massachusetts General Hospital. Initially, we acquired the subject's 
Tl-weighted MRI data (1 mm isotropic) on a Siemens Trio 3 T scanner and 
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reconstructed the tissue compartment boundary surface meshes using the SimNIBS, 
and Freesurfer tools [34, 35]. The subject’s anatomical landmarks were registered to 
the MRI data using the Localite neuronavigation system and two separate trackers 
were mounted on the 3-axis coils to track their position with respect to subject’s 
head. We set up a TCP/IP protocol to communicate the coordinates of these coils 
into MATLAB and calculated the combinatory total E-field of the two 3-axis coils. 
A similar setup was used for the commercial C-B60 TMS coil that was used as a 
reference. 

The EMG data was continuously recorded from the First Dorsolateral interosse- 
ous (FDI) muscle of the right hand using the BrainAmp ExG amplifier (Brain 
Products, Gilching, Germany) and the Motor Evoked Potentials (MEPs) corre- 
sponding to each TMS pulse were saved for post-hoc analysis. Either one or two 
commercially available TMS stimulators were used (MagPro X100, MagVenture) 
to deliver the current pulses. To identify the FDI MEP hotspot, we started with the 
stimulation intensity set to 50% MSO and moved the C-B60 TMS coil over the M1 
cortex guided with the neuronavigation system with the concurrent display of the 
E-field similar to the setup described in Fig. 4. Single pulses were delivered, and the 
intensity was increased until an MEP amplitude of ~50—100 „V was observed. The 
smallest stimulator output intensity by which the EMG response of at least 50 uV 
with 50% probability is observed, was then defined as the FDI resting motor thresh- 
old (MT) [36]. For quantitative reference, we placed the C-B60 coil at the previ- 
ously defined FDI ‘hot spot’ and recorded the MT as well as several clearly 
suprathreshold MEPs (at 107% resting MT). Next, the 2x3-axis coil array was posi- 
tioned over the FDI hot spot and the z-elements of the two 3-axis coils were syn- 
chronously activated to mimic a figure-of-eight coil configuration. The stimulation 
intensity was increased until clearly suprathreshold MEPs were observed. The post- 
hoc analysis of the MEP responses in Fig. 5 shows that both dual z-elements and the 
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Fig. 5 (a) The z-elements of two 3-axis coils stimulating the left M1 hand area. (b) The recorded 
MEPs from FDI muscle. (c) The E-field distribution on the white matter surface. (d, e) 
Corresponding results for the C-B60 coil 
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C-B60 coil elicited suprathreshold responses with similar EMG morphologies 
recorded at the FDI muscle. Furthermore, the computationally estimated cortical 
E-fields show robust agreement for both the amplitudes and the spatial patterns. 


6 Summary and Discussion 


A key requirement of multi-channel TMS technology is fast and accurate computa- 
tion of the E-fields from all coil elements to determine the optimal linear combina- 
tion of the input current values in order to stimulate a desired cortical target region. 
In this study, we demonstrated that precise and fast computation of the E-field can 
be achieved by (1) reconstructing high resolution surfaces meshes based on the 
individual MRI data, (2) reducing the computational burden of high resolution 
E-field modeling by leveraging the MSP-based approach utilizing a spatially fixed 
dipole basis set, and (3) retrieving the pre-calculated solutions for the near-real time 
step to obtain the total E-field by means of simple matrix multiplications based on 
matching the incident field of the TMS coil with the dipole basis set [19]. The posi- 
tion of the TMS coils with respect to subject’s head can be acquired with an optical 
tracking camera and displayed in a neuronavigation software with frame rates of 
10-15 Hz. This gives us a time frame of 60-100 ms for E-field calculations to be 
performed for an interactive navigation pipeline. Our results suggest that we are 
able to calculate the E-field of a C-B60 coil within this timeframe on a high- 
resolution WM surface with about 100 K triangular mesh elements. For the multi- 
channel TMS arrays the E-field computation time will increase based on the number 
of coils and coil elements used. However, we plan to improve the computational 
efficiency further by utilizing a parallel computation approach to obtain the E-field 
of each coil in the array. 

Additionally, we showed that the z-elements of the two 3-axis coils are capable 
of clearly suprathreshold stimulation when used in a ‘dual-channel array’ driven by 
synchronized triggering of two independent TMS stimulators. For more flexible 
stimulation, the orthogonal placement of the elements in the 3-axis coil provides 
efficient decoupling between the elements while granting three degrees of freedom 
in shaping the induced E-field pattern. Moreover, arranging several 3-axis coils in 
an array will further increase the degrees of freedom for efficient stimulation of a 
target area with a desired field orientation without the need for physical movement 
of the coils. We have recently developed a novel 9-channel stimulator system in 
partnership by MagVenture, by which each element of the 3-axis coils can be driven 
by independent stimulator units, allowing delivery of pulses with high com- 
bined power. 

The experimental results show that the juxtaposition of two Z-elements provides 
an E-field intensity and spatial distribution similar to the commercially available 
C-B60 coil, producing clear suprathreshold stimulation as measured with MEPs. 
Furthermore, the MSP-approach used in multi-channel TMS array can be used as a 
pre-planning step to rapidly calculate the E-field at several locations based on 
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specific criteria of the experiment [37] or to optimize the positions/orientations of 
the individual 3-axis coil elements. The biophysical signals recorded during the 
experiment can also be paired with the previously calculated E-field patterns for 
post-hoc analysis [38, 39], or be utilized as feedbacks to adjust the stimulation 
parameters (e.g., input current to each coil) in a closed-loop stimulation paradigm 
[40—42]. 
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1 Introduction 


Magnetic and electrical brain stimulation therapies are widely used to treat neuro- 
degenerative disorders. One of the commonly used non-invasive techniques is tran- 
scranial magnetic stimulation (TMS) that employs magnetic induction to stimulate 
the brain to improve symptoms for diseases such as depression. Computational 
modeling has been used to assess the effectiveness and safety of TMS. A detailed 
brain model is available to allow for these assessments [5]; however, the model is 
only based on one subject. In order to allow for careful planning of a given treat- 
ment regimen, modeling needs to be completed on a per-patient basis. 

A patient-specific brain model can be created using medical imaging data. 
Magnetic Resonance Imaging (MRI) structural data is often used for such purpose. 
There are a variety of semi-automatic segmentation methods available [1, 2, 4, 22] 
that can generate a 3D head model using a set of T,- and T -weighted images. 
Segmentation varies across different methods [6, 18]. Therefore, it may affect elec- 
tric field distributions in electric-field modeling. Some of the segmentation methods 
resulted in differences in electric field distributions of up to 30% when evaluated 
with one computational modeling method [11]. More investigations are needed to 
confirm the degree of differences among different image segmentation and compu- 
tational modeling methods. 

In this study, the T,- and T»-weighted images of 16 subjects were processed with 
three different segmentation methods. Computational modeling of TMS was per- 
formed based on each segmented data by targeting both the primary motor cortex 
and the dorsolateral left prefrontal cortex (DLPFC), then the simulated electric field 
results were compared and evaluated. 
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2 Materials and Methods 


2.1 MRI Data and Segmentation 


MRI T,- and 7>-weighted images were used from 16 Human Connectome Project 
healthy subjects [19] with an isotropic resolution of 0.7 mm per voxel. Two pipe- 
lines implemented in the SimNIBS software package v3.2 [15] were used for seg- 
mentation, headreco [10] and mri2mesh [20], as well as a highresolution FreeSurfer 
[2] pipeline (fshires) [21]. Both the headreco and mri2mesh segmentation methods 
generate surface and volume segmentation of brain and head structures including 
gray matter (GM), white matter (WM), cerebrospinal fluid (CSF), skull, and skin. 
The fshires pipeline produces GM and WM segmentation based on the native sub- 
millimeter resolution. 

The surface meshes were generated using the default options from SimNIBs for 
both headreco and mri2mesh and the high resolution option within Freesurfer 
(fshires) using the -hires flag. The default surface resolution yields surface meshes 
containing a combined total of roughly 800,000 to 1 million facets. 


2.2 Electromagnetic Simulation 


A boundary element fast multipole method (BEM-FMM ) solver was used for elec- 
tromagnetic modeling [7—9]. The solver utilizes the generated surface meshes for 
field estimation. A figure-eight TMS coil was modeled with a diameter of 90 mm for 
each loop. The coil model was modeled based on a commercial coil (MRiB91 of 
MagVenture, Denmark). The coil was placed to target both the patient's left primary 
motor cortex (the hand knob) and the DLPFC via a projection approach and sulcus- 
aligned mapping [3, 12]. These regions were chosen because they are common tar- 
gets for TMS therapy. An example positioning of the coil is shown in Fig. 1. Two 
target points were used for each subject, located within the primary motor cortex 
and the DLPFC, respectively. The coil position was determined using three steps. 
First, the coil was placed so that the centerline (shown as the black line in Fig. 1) 
passed through the given target point on the gray matter interface. Second, the coil 
centerline was made to be perpendicular to the skin surface. Lastly, the coil position 
was adjusted so that the dominant field direction was roughly perpendicular to the 
nearest sulci [7]. 


2.3 Analysis 


Average surface displacement between mesh surfaces generated based on the head- 
reco and mri2mesh segmentation methods were compared across the 16 subjects for 
the gray matter, white matter and CSF surfaces. The displacement was calculated by 
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Fig. 1 An example TMS Coil with Gray Matter 
placement of the TMS coil 


targeting the subject’s left 
motor cortex is shown. The 
black line (the coil axis) 
was used to confirm the 
coil placement and runs 
perpendicular to the coil 


taking the mean of the shortest distance from every triangle centroid of the relevant 
surface from one segmentation method to all triangle centroids of the surface from 
the other segmentation method. The electric field values were compared by extract- 
ing the values in a 100 mm line perpendicular to the TMS coil axis, along the black 
line shown in Fig. 1. Comparisons were performed in pairs between mri2mesh and 
headreco and between mri2mesh and fshires. These values were extracted for both 
of the target points for each subject. The average electric field difference, maximum 
absolute difference, and maximum percentage difference were compared for all 16 
subjects. The average field differences at the target point for both the motor cortex 
and DLPFC were tested for statistical significance using a paired t-test. 

Finally, for additional visualization of the electric field stimulation mapping, all 
subject electric field results were mapped onto the inflated surface from the FreeSurfer 
common space (fsaverage). The electric field was exported and mapped onto this 
common surface to map average electric field and its difference surface over the 16 
subjects. These maps allow for the qualitative evaluation of focality variations 
between three different segmentation methods over the entire subject space. 


3 Results 


All extracerebral and cortical surfaces were successfully reconstructed using the 
three segmentation methods. The average surface displacement was only calculated 
for regions of the brain located in the superior cerebral cortex for all 16 subjects 
between headreco and mri2mesh and the results are summarized in Fig. 2. The dis- 
placement between fshires and mri2mesh were not compared as the surfaces use the 
same algorithm from FreeSurfer and the comparison was done by other study [21]. 
The CSF surface displacement was three times more than the white matter surface 
displacement, with an average difference of 0.9 mm (+0.2 mm) for the CSF and 
0.3 mm (+0.06 mm) for the white matter. An example of all the surfaces overlaid on 
a subject T;-weighted image is shown in Fig. 3. 
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Fig. 2 Average displacement between headreco and mri2mesh surfaces in millimeters across all 
16 subjects for CSF, GM, and WM 


For the target point located in the motor cortex, the electric fields showed similar 
distributions for the mri2mesh and the headreco segmentation methods. There was 
an average difference in magnitude of 0.8 V/m and a maximum difference of 
61 V/m. The average percentage difference was 2%. In the region of interest within 
5 mm of the target point, the average percentage difference was approximately 596. 
For the mri2mesh and fshires segmentation methods, the average percentage differ- 
ence in the region of interest surrounding the target point was 0.796. Extracted elec- 
tric fields (along the dotted line) and the surface contour lines are shown on the 
subject's T, image in Fig. 3. 

For the target point located in the DLFPC, there was an average difference of 
0.8 V/m in magnitude and a maximum difference of 50 V/m between the mri2mesh 
and headreco segmentation methods. The average percentage difference was 2.8%. 
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Fig. 3 (continued) field results and the thin contours show the headreco field results (fshires 
contours are not shown because they are almost identical to the mri2mesh results). Color indica- 
tions for surfaces are red for skin, orange for skull, yellow for CSF, cyan for gray matter, and 
purple for white matter. The dotted white line on the axial cross section is a projection onto the XY 
plane of the 100 mm line running along the axis perpendicular from the coil. The electric field 
result along the dotted white line was extracted as shown at the bottom of the figure. The dotted 
black line indicates the location of the target point of stimulation. The field peaks were observed at 
the anatomical structure transition points (represented by arrows): (1) Skin-Skull, (2) Skull-CSF, 
(3) CSF-GM, (4) GM-WM 
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Fig.3 Example of electric field targeted on motor cortex is shown along with an axial cross sec- 
tion of a subject (target stimulation point in magenta). The thick contours show the mri2mesh field 
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In the region of interest within 5 mm of the target point, the average percentage dif- 
ference was approximately 1.8%. For the mri2mesh and fshires segmentation com- 
parison, the average percentage difference in the region of interest surrounding the 
target point was 0.1%. Extracted electric fields (along the dotted line) and the sur- 
face contour lines are shown on the subject’s T, image in Fig. 4. 

The electric field difference across all 16 subjects at the target point between 
headreco and mri2mesh was statistically significant (p = 0.005) for the motor cortex 
and not significant between mri2mesh and fshires (p = 0.83). The electric field dif- 
ference at the target point for the DLFPC was not significant for either the compari- 
son between mri2mesh and headreco (p = 0.19) or between mri2mesh and the 
fshires (p = 0.23). The average percentage differences in the electric field over the 
100 mm line for all subjects and target points are shown in Tables 1 and 2. 

The average surface mappings for both the frontal and motor cortices are shown 
in Figs. 5 and 6. The average difference between the headreco and mri2mesh elec- 
tric field results for the motor cortex are shown in Fig. 7. The electric fields and 
corresponding differences are mapped onto the fsaverage inflated surface along 
with cortical parcellation contours. 


4 Discussion and Conclusion 


This study focused on an evaluation of the electric field differences between The 
surface displacement between mri2mesh and headreco were seen at the CSF bound- 
ary, where the CSF surface was estimated closer to the gray matter for the headreco 
segmentation. Such surface displacement aligns with results shown in previous 
studies (cf., [11, 13]). In particular, Seiger et al. 2018 demonstrated that Freesurfer 
was more accurate in its calculation of cortical thickness; however, CAT12 based 
methods (such as headreco) were faster and yielded reliable results [16]. 

Differences between the mri2mesh and fshires were subtle as the underlying 
algorithm within FreeSurfer to segment the T; image is same for both methods. The 
lack of significant differences seen in the electric. 

field between the two methods indicates that additional proceces to use a native 
submillimeter resolution is not necessary for BEM-FMM based computational 
modeling. Nevertheless, this may not be the case when highly detailed submillime- 
ter surfaces are needed such as a very small structure. 

In both the motor cortex and the DLPFC, the low average percent difference in 
the electric field suggests that the effect of the segmentation method differences was 


T ————————————————————————————————————————————— ws 
Fig. 4 (continued) because they are almost identical to the mri2mesh results). Color indications for 
surfaces are red for skin, orange for skull, yellow for CSF, cyan for gray matter, and purple for 
white matter. The dotted white line on the axial cross section is a projection onto the XY plane of 
the 100 mm line running along the axis perpendicular from the coil. The electric field result along 
the dotted white line was extracted as shown at the bottom of the figure. The dotted black line 
indicates the location of the target point of stimulation. The field peaks were observed at the ana- 
tomical structure transition points (represented by arrows): (1) Skin-Skull, (2) Skull-CSF, (3) 
CSF-GM, (4) GM-WM 
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Fig. 4 Example of electric field targeted on DLPFC cortex is shown along with an axial cross 


section of a subject (target stimulation point in magenta). The thick contours show the mri2mesh 
field results and the thin contours show the headreco field results (fshires contours are not shown 
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Fig. 5 Mean electric field values are mapped on the inflated surface for the motor cortex stimula- 
tion across all 16 subjects for mri2mesh (a), fshires (b), and headreco (c) along with the cortical 
parcellation contours. All three methods showed high electric field values in the target (precentral 
gyrus) along with the postcentral and caudal middle frontal gyri. No notable differences were 
Observed between fshires and mri2mesh, and higher electric field values were observed in all three 
regions for the headreco method 


35 Vim 70 Vim 


Fig. 6 Mean electric field values are mapped on the inflated surface for the DLPFC stimulation 
across all 16 subjects for mri2mesh (a), fshires (b), and headreco (c) along with the cortical parcel- 
lation contours. All three methods showed high electric field values in the target (superior frontal 
gyrus). No notable differences were observed between fshires and mri2mesh, and slightly lower 
electric fields were observed in all three regions for the headreco method 


minimal. Although the overall difference was low (j 5%), the localized field differ- 
ence near the target point across all 16 subjects was statistically significant for the 
motor cortex and could affect the intended stimulation there. The field difference for 
the DLPFC was not significant and the fields were more similar between headreco 
and mri2mesh as shown in Fig. 4. 
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Motor Cortex Difference Mapping 


Fig. 7 Mean electric field differences between mri2mesh and headreco on the inflated surface are 
shown for the motor cortex stimulation. The area where the headreco electric field was lower than 
the mri2mesh electric field is shown in blue, and the area where the headreco electric field was 
higher than the mri2mesh electric field is shown in red 


The percent difference in the electric field for the DLPFC was similar to that of 
the motor cortex along the entire 100 mm line for all 16 subjects. However, the 
region directly surrounding the target point showed lower variability for the DLPFC 
target compared to the motor cortex (Fig. 4). This trend is consistent in the average 
percent difference in the 5 mm region surrounding the target point. Therefore, there 
may be more segmentation variability in the motor cortex. 

Analysis of the extracted electric field in Figs. 3 and 4 showed sudden changes 
in the field that resulted from the anatomical structure transitions. Some peaks are 
not aligned at the distance from the coil. For example, the Skull-CSF transitions in 
both figures were approximately 23 mm for the headreco segmentation whereas it 
was approximately 22.5 mm from the coil for the simulation with the mri2mesh and 
fshires methods. These sudden electric field changes resulting from segmentation 
differences can affect TMS therapy because the resulting neuronal excitation is a 
function of the electric field gradient rather than the electric field magnitude. 
Additional work that directly evaluates the gradient of the electric field will provide 
more insight into the effect of the segmentation on the TMS therapy. 

The average electric field mapped on the fsaverage surface showed minimal dif- 
ferences between the fshires and mri2mesh segmentation methods as the differences 
were within 2 V/m. The electric field difference mapped on the fsaverage surface 
revealed that there were clusters that exceed 10 V/m of the electric field differences 
between headreco and mri2mesh (Fig. 7). The differences were also largely one- 
sided, with the fields from the headreco segmentation consistently higher than those 
from mri2mesh. 
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Coil positioning was critical to the electric field estimation. In this study, coil 
positioning was determined automatically based on the topology of the input meshes 
used. The pre-processing steps to select the proper coil position may differ signifi- 
cantly between segmentation methods. 

The resulting field differences for each coil position were small but measurable. 
Though average percent differences observed along the 100 mm observation line 
were less than 5%, significant differences in the electric fields between segmenta- 
tion methods were observed for the motor cortex simulation. Moreover, the field 
differences shown by subtle segmentation differences indicate an importance of 
patientspecific modeling as various previous studies have shown morphometric dif- 
ferences across age [14] and sex [17]. Future studies with different types of TMS 
coils and different segmentation and computational modeling methods may further 
improve a modeling approach for robust treatment for TMS and other neuromodula- 
tion devices. 
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1 Introduction 


Transcranial magnetic stimulation (TMS) has been approved by the Food and Drug 
Administration (FDA) for treatment-resistant major depression [1] and Obsessive- 
Compulsive Disorder [2]. Its therapeutic effects in other psychiatric and neurologi- 
cal disorders, including drug addiction, are emerging [3, 4]. From both clinical and 
basic neuroscience perspectives, there has been a strong demand for stimulation 
tools that can reach deep brain regions with small size targeted stimulations. For 
example, decades of neuroimaging studies have identified malfunction of dorsal 
anterior cingulate cortex, insular and amygdala in a range of psychiatric disorders. 
These structures are 4 cm or more below the scalp. Unfortunately, with current 
technologies, the stimulation targets are limited to superficial brain regions, or oth- 
erwise wide brain areas are stimulated when a deep brain structure is targeted. 
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The output of a TMS coil can be treated as field emission from a finite size aper- 
ture and follows a specific depth-focality tradeoff rule. Deng et al. theoretically 
calculated the depth-focality profiles of 50 TMS coils [5]. Two groups of mainstream 
coils, circular and the figure-8, formed two depth-focality tradeoff curves, respec- 
tively. The study concluded that at shorter depth, which is smaller than 3.5 cm, all 
figure-8 type of coils follows a better depth-focality tradeoff rule and it will be advan- 
tageous to use the figure-8 coil. A number of studies have attempted to design TMS 
coils for enhanced the penetration depth or improved the focality. Rastogi modified 
conventional figure-8 coil to improve its focality but that, on the other hand, signifi- 
cantly weakened the electric field strength generated in brain tissues [6]. Crowther 
suggested a “Halo coil" design to improve the penetration depth of the conventional 
circular coils [7, 8], but this design sacrifices the coil's focality. Luiz modeled multi- 
channel coil arrays to improve the focality and penetration depth profile [9]. However, 
this design involved complicated coil structures, and it required higher efficiency of 
the coils’ cooling system. Alternative coil design strategy is needed to go beyond the 
depth-focality tradeoff limitation. Roth et al. have developed the H-coil for human 
deep brain stimulation, but the design is still limited by the depth-focality tradeoff 
with a relatively large field spread [5, 10, 11]. We recently reported a multi-layer 
winding-tilted coil design approach for focused rodents’ brain stimulation that 
induced unilateral movements [12]. The goal of this study is to extend this novel 
strategy to design TMS coils for deep and focused human brain stimulation and 
compare its depth-focality characteristic with other conventional coils. 


2 Methods 


Figure la illustrates our TMS coil design (see reference [12] for details). In this 
study, the coil dimensions have been extended for human brain stimulation. For 
each single circular coil, the inner and outer diameters are 8 cm and 9 cm, respec- 
tively. The thickness of each single coil is 1 cm and 5 circular coils are accumulated 
along the central axis in this model. Considering the overall inductance of the coil 
model, we further extended the coil length while using narrow (lower value of “d,— 
d;’) but thick (higher value of h) coil windings, so that the total turn number could 
be low enough to control the coil inductance. 

The head model is a homogenous sphere with the diameter of 17 cm and isotro- 
pic electrical conductivity of 0.33 S/m^!, The definitions of stimulation depth and 
focality in Fig. 1b are based on the half-value depth (di2) and half-value volume 
(Vi). Modeling frequency was 5 kHz All conditions are identical to the modeling 
in the study by Deng et al. [5], except that the software we used for calculation is the 
COMSOL AC/DC module (finite element analysis software, COMSOL Inc.), which 
is different from that in their study. To calibrate our calculation using COMSOL, we 
firstly selected 4 coils, which were already documented in the study by Deng et al., 
and compared our results with theirs. The selected coil models for calculation cali- 
bration were the 50 mm, 70 mm and 90 mm circular and the 70 mm figure-8 
Magstim coils. 
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Fig. 1 (a) Simulation model and induced electric field distribution in the human head model; (b) 
The definitions of the stimulation depth and focality 


Table 1 Comparison of the stimulation depth (d;;) and focality (Sj) calculated by COMSOL 
with data recorded in Deng's study 


50mm circular |70 mm circular |90 mm circular |70 mm figure-8 

Coil type magstim coil magstim coil magstim coil magstim coil 
dy, (by 1.31 1.49 1.75 1.45 
COMSOL) / cm 
dy, (by Deng 1.29 1.44 1.74 1.41 
et al.) / cm 
S 4 (by 53.1 65.8 87.8 13.8 
COMSOL) / cm? 
S y (by Deng 53.7 66.0 87.4 14.8 
et al.) / cm? 

3 Results 


3.1  Cross-Validation of Theoretical Simulation 


The stimulation depth and focality calculated by COMSOL of the 4 selected coils 
(50 mm, 70 mm and 90 mm circular and the 70 mm figure-8 Magstim coils) as cali- 
brations are listed in Table 1, compared with the data in Deng's study [5]. Both the 
stimulation depth and focality from the two finite element analysis software matched 
reasonably well with minor differences. 

Theoretical simulation of multi-layer winding tilted coil design. 
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For our multi-layer winding tilted coil design, we firstly investigated how the tilt 
angle O affected the stimulation depth and focality. The angle was adjusted from 0 
degree to 70 degrees with a step of 10 degrees without changing any other parameters 
in the modeling. The values of the stimulation depth and focality were marked in the 
same plot summarized in Deng et al. in Fig. 2 [5]. Coils with 0-degree tilt angle (or flat 
coils) are located along the circular coil curve. As the tilt angle increases, the location 
of the coil in the depth-focality tradeoff plot moves from the circular coil curve towards 
the figure-8 coil curve. If the number of the winding layer is set to 5, the coil location 
drops on the figure-8 coil curve when the tilt angle reaches 50 degrees. Further enlarg- 
ing the tilt angle enables the curve, which is plotted from the trace of our coil designs 
(formed by the blue square dots in Fig. 2), to penetrate the figure-8 coil curve. For 
example, when the tilt angle is adjusted to 60 or 70 degrees, the locations of the coils 
in the plot are below the figure-8 coil curve, and that indicates a better depth-focality 
characteristics than the existing TMS coils. The number of the winding layers are also 
adjusted to 2 and 9. We finalize that a smaller number of winding layers moves the 
curve of the coil in the plot towards the left side. For a certain tilt angle, both the stimu- 
lation depth and focality decrease. However, when the winding layer number increases 
from 5 to 9, the stimulation depth is not considerably improved. This phenomenon 
may be cause by the longer distance from the stimulation target to the few top layers, 
which significantly weakens the electric field strength at the stimulation target contrib- 
uted by those layers. Figure 3 presents the induced electric field distribution on the 
human head model surface by the proposed coil designs with various design parameters. 


eCalibrated Data x2 Winding Layers m5 Winding Layers 49 Winding Layers 


400 


S'4 (cm?) 
E 


4 
0.4 d% (cm) 2 


Fig.2 Calculations of stimulation depth and focality for the multi-layer winding-tilted coil design 
with air core and tilt angle ranging from 0 to 70 degrees (green curve), the number of winding lay- 
ers of 2, 5 and 9 
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To further demonstrate the advantage of our coil design, we compared the elec- 
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Fig. 3 Induced electric field distribution on the human head model surface by coil designs with 
different design parameters 


tric field decay rate in the brain model for 6 coil structures, as shown in Fig. 4, and 
experimentally demonstrated the tilt angle could improve the focality of the induced 
electric field distribution. The electric field decay rate curves in the human head 
model in Fig. 4 indicate that neither the application of ferromagnetic core nor the tilt 
angle 0 is able to improve the field decay rate, but the winding accumulation along 
the coil's central axis considerably improves it. For example, at the depth of 3 cm, a 
5-layer winding accumulation improves the electric field decay rate by 4-5%. 


3.2 Experimental Validation 


To verify how the angle 0 affected the focality, we fabricated 3 coil prototypes with 
winding's tilt angles of 20, 10 and 0 degree. They shared the same coil length, inner 
and outer diameters, which were 4.4 cm, 3.8 cm and 7.5 cm respectively. Each coil 
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Fig. 4 Comparison of field decay rates for 6 different coil structure models: (a) single-layer flat 
circular coil with air core; (b) 5-layer flat circular coil with air core; (c) single-layer circular coil 
with air core and 40-degree tilt angle; (d) 5-layer circular coil with air core and 40-degree tilt 
angle; (e) single-layer flat circular coil with ferromagnetic core; (f) 5-layer flat circular coil with 
ferromagnetic core 


was wrapped by 20 turns of the litz wires, and each turn contained a bundle of 135 
piece of AWG30 wires. The TMS coil was driven by a customized driving circuit, in 
which an insulated gate bipolar transistor (IGBT) was used as the switch to control 
TMS pulses and a capacitor bank charged by a power supply was the main current 
source [13]. The charging voltage of the capacitor bank was set to 100 V and the 
pulse duration was 250 ps. The induced electric field was measured with a modified 
Rogowski coil electric field probe customized in our lab [14]. The electric field was 
mapped in the medium of air within a plane 2 cm away from the coil surface. Since 
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Fig. 5 Induced electric field measurements using modified Rogowski coil probe for the 3 coil 
prototypes with 20 degrees, 10 degrees and 0 degree wire wrapping tilt angle at the depth of 2 cm 
in air medium 


the electric field components along the Z axis (in parallel with the coil’s central axis) 
was small enough to be negligible, only the X and Y components of the electric field 
were measured. Considering the probe size, we mapped the electric field at a step 
size of 5 mm within the X-Y plane. An area of 8 cm x 8 cm was scanned for each coil. 

Figure 5 shows the heat maps of the electric field distributions for the 3 scanned 
coils. For the coil with a tilt angle of 20 degrees, the area with the field strength over 
or equal to 80% of its peak value (Epeak) is only 14.5 cm’; while for the other 2 coils, 
the values have reached 22.5 cm? (for 10-degree tilt) and 30.75 cm? (for 0-degree 
tilt), respectively. It is also found that a larger tilt angle of the coil windings would 
slightly increase the coil’s inductance. For example, the inductances of the 3 scanned 
coils were measured to be 52.2 uH, 51.21 uH and 49.4 uH, respectively. 


4 Discussions 


The combination of air-core accumulated windings along the coil's central axis and the 
tilt angle of the windings provides a significant innovation to the depth-focality profile 
of TMS coils. The tilt angle technique is a mild symmetry breaking method. It does not 
reduce the equivalent field emission aperture size or speed up the electric field decay 
rate along the coil’s central axis direction, but significantly distorts the ring shape elec- 
tric field distribution, resulting in a much smaller focal spot. The only limitation of this 
coil design, to our knowledge, is its possible larger inductance compared with conven- 
tional TMS coils. The current flowing inside the coil I(t) can be expressed as. 


I(t) =~ sin(«t)exp(-o"), (1) 


where œ and o are related to the charging of the capacitors (C) in the stimulator cir- 
cuit, the inductance of the coil (L), and the resistance in the LC circuit. V, is the volt- 
age to charge the capacitors [15]. The induced electric field E(t) can be expressed as 
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E(r)-a 1 za £. (2) 


So, E is inversely proportional to L. The Magstim human TMS coils always have an 
inductance of around 20 uH (from the Magstim Rapid [2] system manual) [16]. Our 
prototype human TMS coil design may have a higher inductance due to its length. 
This would require a stimulator of higher power output to drive the coil and enhance 
the current load in it, and on the other hand, the turn number of the coil can be 
reduced to limit its inductance to a reasonable value. 

The curves we plotted for our multi-layer winding-tilted coil design in the depth- 
focality tradeoff profile in Fig. 2 demonstrate better focality than the H coils. The 
half-value depth values are comparable with the H coil designs. The half-value 
depth values are larger than the conventional figire-8 coils. Our coil design has 
achieved better depth-focality characteristic than a large sized double cone coil for 
human brain stimulation [5]. However, the double cone coil is known to induce 
scalp and facial pain and has limitations on its stimulation targeting site due to its 
geometry [17]. Our coil design has the advantage of much smaller contact area with 
the human head during the stimulation, and that may reduce or even exempt the 
scalp or facial pain. Moreover, with our coil design, it is feasible for users to conduct 
multisite stimulation, which cannot be accomplished by the double cone coils. Our 
design has provided the current best method for multisite human brain stimulation 
than all the conventional TMS coils considering both the stimulation depth and 
focality. The stimulation depth and focality are adjustable by the geometry of the 
coil design, for example, by tuning the tilt angle of the windings and the coil's outer 
diameter. This novel design provides a promising solution for the future deep and 
focused multisite human brain stimulation. 
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in Transcutaneous Spinal Direct Current 
Stimulation (tsDCS) 
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Mamede de Carvalho, and Pedro C. Miranda 


1 Introduction 


Transcutaneous spinal direct current stimulation (tsDCS) is a non-invasive stimu- 
lation technique considered as a possible therapeutic resource for spinal cord dys- 
functions such as chronic pain or motor system lesion [1]. tsDCS consists in the 
application of low intensity direct currents of 2—4 mA to the spinal cord (SC) using 
electrodes placed over the vertebral column near target regions. Exploratory 
experimental studies in humans pursued in the last decade demonstrated the neu- 
romodulatory potential of tsDCS to change signal transmission along nociceptive 
ascending pathways and spinal motor and reflex circuits in cervical and thoracic 
spinal segments [2-5]. 
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The effects of tsDCS rely mainly on the electric field (EF) induced in the nervous 
tissue, similarly to brain stimulation techniques. The EFs may contribute to inhibit 
or facilitate neuronal responses, by transiently changing the resting membrane 
potential. This effect will depend on how each spinal neuron or circuit is orientated 
relative to the EF. Just as in brain stimulation, the spatial distribution of the EF 
induced by tsDCS depends on electrode number, design (shape and structure), 
placement relative to target, current intensity and polarity (anodal/cathodal). 
Computational studies using numerical methods to predict the EF and current dis- 
tribution in realistic human models may be powerful tools to fine-tune experimental 
tsDCS protocols targeting specific spinal regions of interest [6—9]. 

Modelling studies in tsDCS generally consider square or rectangular electrodes 
in bipolar montages, which originate EFs in the spinal cord that increase with dis- 
tance between electrodes at the cost of a larger (less focal) target region. Cortical 
stimulation studies predict that montages using smaller electrodes in a multi- 
electrode setup can induce maximum EF in a smaller target region [10-12]. Thus, 
one question of interest is to determine if a multi-small-electrodes paradigm in 
tsDCS can also originate a more focal EF in the spinal cord, as in brain stimulation, 
which could be relevant for neuromodulation of specific spinal circuits. 

The effects of different montage settings are also dependent of the electrical 
properties of tissues located in the current path to the stimulation target. Sensitivity 
analyses performed on transcranial direct current stimulation (tDCS) demonstrated 
changes of the EF magnitude when varying the conductivities of tissues. Scalp, 
CSF, gray matter (GM) and skull were the tissues that introduced larger variability 
on EF values in the cortex [13, 14]. The EF induced in tsDCS also show spatial 
characteristics, such as local hotspots, that seem to be related with anatomical fea- 
tures such as vertebrae edges and disks protrusions into the spinal canal. These local 
features may be related with an interplay between the conductivities of different 
tissues, such as vertebrae, intervertebral disks and CSF [7-9]. 

This chapter is dedicated to study the interplay between the conductivities of tis- 
sues and two different types of electrode montages on the spatial distribution of the 
EF induced by tsDCS. These types of montages will be modelled considering tsDCS 
over the lower thoraco-lumbar spinal cord: a montage with two large and square 
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electrodes (50x50 mm?) and a grid with 4x4 small circular electrodes (@ = 10 mm). 
For these two types of montages, EF changes will be investigated considering rota- 
tion of the grid or increases in distance between the two large electrodes. The impact 
of electrical conductivities of the tissues located between the electrodes and the SC 
will be addressed in each montage considering three different types of models: 
homogeneous, semi-homogeneous and heterogeneous. The clinical relevance of 
differences in the EF induced by the multi-electrode grid system and two-electrode 
montages in tsDCS application will be discussed. Furthermore, the influence of dif- 
ferent conductivity properties of tissues will be studied since it can be determinant 
to optimize stimulation delivery for intended spinal targets. 


2 Methods 


2.1 Realistic Human Model and Electrode Placement 


This study used the same realistic human model as in previous works (e.g [9]). This 
model was based on relevant tissue masks from the 34 years-old male Duke model 
of the Virtual Population Family [15], comprising 15 tissues (Fig. 1). The spinal GM 
was artificially designed considering general anatomical knowledge and measure- 
ments from the Visible Human Data Set (National Library of Medicine, NLM, 
Visible Human Project®, www.nlm.nih.gov/research/visible/visible human.html). 
The full model was truncated at the level of the thighs and above the elbows, to 
reduce model size and computational time. 

Electrodes were represented as hydrogel layers with 1 mm thickness, consider- 
ing two geometries available for standard electrodes used in transcutaneous electri- 
cal nerve stimulation (TENS): 50 x 50 mm? square electrodes (Fig. 2a); circular 
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c(skin) = 0.435 S/m o(fat)=0.040S/m — c(muscie) = 0.355 S/m (av) a(bone) = 0.006 S/m aivertebrae) = 0.006 S/m — c(cerebellum) =0.290S/m erat 
e(lungs) = 0.046 S/m (av) o(disks) = 0.200 S/m (brainstem) = 0.154 S/m 
o(heart) = 0.535 S/m (av) ^ c(dura) = 0.030 S/m G(spinal-WM) = 0.143 S/m 
olviscera) = 0.123 S/m (av) — c(CSF) = 1.790 S/m o(spinal-GM) = 0.333 S/m 


Fig. 1 Tissues included in the human model and corresponding isotropic values of electrical con- 
ductivity c, as assessed in Fernandes et al. [9]. (av) indicates that conductivities were averaged 
over different tissues' components or longitudinal/transversal values (CSF — cerebrospinal fluid; 
WM - white matter; GM — gray matter) 
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Fig. 2 Geometry and dimensions of the square (a) and circular (b) electrodes represented in the 
model. Variations considered for 2-electrode montages (c) and Grid montages (d); anodes and 
cathodes are represented by red and blue electrodes, respectively 


electrodes with a diameter of 10 mm (Fig. 2b). TENS electrode design was consid- 
ered in this study to adequately reproduce the electrodes that will compose the grid 
montages to be studied by one of our clinical teams. Two types of electrode settings 
were considered: bipolar montages with two square electrodes (2-electrode); multi- 
electrode montages with 16 circular electrodes placed in a 4x4 electrode array 
(Grid). Both montages were centered midway between T9 and T11 vertebrae spi- 
nous processes. The 2-electrode montage was varied in alignment relative to the 
vertebral column — transversal (2-T), diagonal (2-D), longitudinal (2-L) — and dis- 
tance — 5 mm and 65 mm between electrode edges (Fig. 2c). The electrode Grid was 
rotated by 45? and 90? relative to the first position (rotation 0°, Fig. 2d). Surface 
meshes were optimized and assembled, and volume meshing was performed with 
the 3-MATIC module from MIMICS (MIMICS software, v16), resulting in 2x10" 
tetrahedral elements for the entire human model with electrodes. Average meshing 
time was 4 hours per electrode montage. 
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Table 1 Isotropic electrical conductivity models 


Model Name c values 
Homogeneous Homogeneous 6 — 0.289 S/m for all tissues 
Semi-homogeneous GM o(GM) = 0.333 S/m; o(other 
tissues) = 0.289 S/m 
WM o(WM) = 0.143 S/m; o(other 
tissues) = 0.289 S/m 
CSF o(CSF) = 1.79 S/m; o(other 
tissues) = 0.289 S/m 
Dura o(dura) = 0.030 S/m; o(other 
tissues) = 0.289 S/m 
Vertebrae o(vertebrae) = 0.006 S/m; o(other 
tissues) = 0.289 S/m 
Disks o(disks) = 0.200 S/m; o(other 
tissues) = 0.289 S/m 
Heterogeneous Isotropic c values from Fig. 1 


GM gray matter, WM white matter, CSF cerebrospinal fluid 


2.2 Electrical Properties of Tissues and Electrodes 


Tissues and hydrogel were assumed to be purely resistive with isotropic electrical 
conductivities. An electrical conductivity of 0.02 S/m was assigned to hydrogel in 
all simulations [16]. Three different levels of tissue heterogeneity were considered 
(Table 1): 


Heterogeneous model — Tissues and hydrogel were assumed to be purely resistive 
with isotropic electrical conductivities, considering DC electrical tissue proper- 
ties compiled in our previous work [9] (Fig. 1); 

Homogeneous model — All tissues were assigned an isotropic electrical conductivity 
of 0.289 S/m, that corresponds to the volume-weighted average of the isotropic 
conductivities of each tissue included in the model; 

Semi-homogeneous model — Considers the same conductivity as the homogeneous 
model for all tissues, except for one tissue in the vertebral column (spinal gray 
matter (GM), spinal white matter (WM), CSF, dura, vertebrae or disks) resulting 
in six different models. 


2.3 Electric Field Calculations 


Electric field (EF) spatial distribution was simulated with COMSOL Multiphysics 
using the finite element method (FEM). Total current intensity was considered as 
4 mA (the maximum delivered by a typical tsDCS device): 4 mA/electrode in 
2-electrode montages and 0.5 mA/electrode in Grid montages. Boundary conditions 
were applied according to [17], considering the hydrogel top surface as isopotential. 
A total of 72 simulations (9 montages variations x 8 conductivity models) were 
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performed, with 2.2x107 degrees of freedom and solution time of about 30 minutes 
per simulation on a computer with 2 quad-core Intel&9 Xeon( processors clocked at 
3.2 GHz and 48 GB of RAM. 

Modelling studies in tDCS reported volume-average E-field values larger than 
0.15 V/m over the hand knob, when reproducing clinical settings with observed 
neuromodulatory effects [17, 18]. tsDCS neuromodulation will be assumed if the 
average EF exceeds this value in the SC. 

EF components were defined as 3 orthogonal vectors: Ejj,, — caudal-rostral ori- 
ented and tangent to SC axis; E44 — ventral-dorsal oriented and perpendicular to SC 
axis; E, — right-left oriented and perpendicular to SC axis. Spatial profiles of the 
magnitude of the total EF and of its components along the spinal WM and GM were 
determined by considering EF averages calculated at 1-mm thick axial slices along 
the SC length. 


3 Results 


3.1 Current Delivery for 2-Electrode and Grid Placements: 
Safety Considerations 


The electrode montages considered in this study may present some safety concerns 
in terms of tissue damage that should be addressed: (1) in 2-electrode placements, 
an inter-electrode distance of 5 mm may be too small to prevent current local max- 
ima at electrode edges; (2) small circular electrodes in Grid montage may deliver 
large current densities below each electrode. 

Current density was predicted at skin and target tissues (spinal GM and WM) for 
the heterogeneous model considering all placements addressed and is summarized 
in Fig. 3. Current density magnitude shows hotspots in skin regions near electrode 
edges (Fig. 3a). For inter-electrode distances of 5 mm, these hotspots are higher 
between electrodes. In Grid montages, the current density is larger at electrodes’ 
edges and maximum in skin regions located below the grid where electrodes’ polar- 
ity changes from anodal to cathodal. Maximum values of the current density vary 
from 6.89 to 24.60 A/m? in skin regions near and below the electrodes, with the 
highest values predicted in Grid montages and the lowest values in 2-electrode mon- 
tages with larger inter-electrode distances (Fig. 3b). The opposite occurs at the tar- 
get tissue, where current density maxima are larger for 2-electrode montages with 
more distance between anode and cathode. 

Current density threshold for cutaneous and nerve lesion was established to be 
143 A/m? for DC stimulation [19]. This value has not been reached in previous 
tsDCS and tDCS experimental studies with no report of tissue damage [1, 5, 20, 21]. 
These studies applied currents of 2-4 mA using larger electrodes at larger distances. 
Models of these tsDCS protocols predicted current densities maxima of 12.6-21.7 A/ 
m? and 0.11—0.15 A/m? at skin and spinal-GM, respectively [9]. These values are of 
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Fig. 3 (a) Current density magnitude distribution in skin below and near electrodes for each mon- 
tage. A color scale is presented at the top right: all distributions are normalized to the respective 
maximum value. (b) Maximum values predicted for the current density magnitude in skin and 
spinal WM and GM for each montage 


the same order of magnitude of the values predicted for the 2-electrode and Grid 
montages and variations presented in this study, thus the montages proposed here 
can be considered safe for future experimental studies. 


3.2 EF Distribution for Grid and 2-Electrode Montages 
in the Heterogeneous Model 


The heterogeneous model includes isotropic electrical conductivities for each tissue 
in the model (Fig. 1), thus providing a more realistic insight on the spatial charac- 
teristics of the EF at target (spinal cord) and surrounding tissues. Figure 4 shows the 
EF magnitude in a target volume of the spinal cord, comprising lower thoracic to 
lumbosacral spinal regions, from T7 to Co segments for all montages and corre- 
sponding variations in rotation, alignment and distance between electrodes. 
2-electrode and Grid montages can be grouped regarding similar EF spatial profiles: 
transversal/rotation 0° (TO); diagonal/rotation 45? (D45); longitudinal/rotation 90° 
(L90). The EF magnitude distribution along the SC is wider and reaches higher EF 
maximum in 2-electrode montages, especially for a larger electrode distance, 
whereas the distribution due to 2-electrode at 5 mm is more similar to that of Grid 
montages (Fig. 4). The EF direction varies with montage groups, which can be seen 
by comparing the maximum values of the EF components, normalized to the maxi- 
mum EF magnitude (Emax; Fig. 4, table at top right). In TO montages, the EF is 
mostly transversal, with a larger right-left component of 0.96-0.98 Emax- In D45 and 
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Fig.4 (a) EF distribution in a selected volume of the SC, between T7 and Co spinal segments for 
all montages and variations; maximum value and name of montage are on top of each plot; location 
of the selected volume is marked by a red rectangle over a representation of the vertebral column 
and colour scale for the EF magnitude are at top and middle left, respectively. (b) Maximum values 
of the EF components normalized to Emax in the spinal WM (top) and GM (bottom). Values of the 
largest component are in blue 


L90, the EF direction is longitudinal, (Ei, (max) = 0.94-0.99 Emax). A maximum 
contribution of E,,,, is present in L90 and 2-D at 65 mm. Maximum values occur 
mainly at T12 and L1 spinal segments, which are located midway between the elec- 
trodes or grid upper and lower borders. Figure 5 shows contributions of different 
montage variations inside the SC, considering an axial slice of the spinal cord at 
T12 level. EF similarities also translate in LO, D45 and L90 grouping at local level. 
D45 and L90 present different spatial patterns in the WM and GM when compared 
to TO. EF orientation projected in the axial plane is also represented in Fig. 5 (black 
arrows) and changes from group to group. D45 presents a EF axial projection with 
a dorsal-ventral, right-left direction, whereas L90 projection is almost only dorsal- 
ventral. This can also be seen when comparing values of E,, and E, normalized to 
maximum in Fig. 4b. This orientation is consistent with anode-cathode relative posi- 
tion in each case. 
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vim 


Fig. 5 EF distribution in an axial slice at T12 spinal level in spinal WM and GM for all montages. 
Left: location of the slice is indicated by a red dashed line over a representation of the vertebral 
column. Right: colour scale, with the red colour corresponding to the maximum EF in the WM in 
each case and orientation of the slices at the top. The black arrows have size proportional to the EF 
magnitude and represent the orientation of the field projected in the axial plane 


Table 2. Maximum values of the EF in the spinal WM and GM in the homogeneous model, in V/m 


2-T 2-D 2-L Grid rotation 
Electrode montage |5 mm |65mm |5mm |65mm |5mm /|65 mm | 0° 45? | 90° 
Spinal-WM 0.62 | 0.69 0.72 | 0.82 0.63 | 0.91 0.51 |0.49 | 0.48 
Spinal-GM 0.56 | 0.64 0.67 | 0.77 0.58 | 0.85 0.46 | 0.43 | 0.43 


3.3 Homogeneous Model: How the Relative Positions 
of Electrodes Influence the EF Direction 


The homogeneous model provides insight on the effect of electrode positions on the 
EF direction and magnitude, reducing variability that arises due to the different 
properties of tissues that surround the SC. 2-electrode montages originate larger EF 
magnitudes in the spinal WM and GM when compared to the 4x4 electrode grid. 
Furthermore, larger distance between the electrode edges results in an EF increase, 
whereas grid rotation does not change the EF considerably in terms of maximum 
values (Table 2). Again, 2-electrode and grid montages can be grouped regarding 
similar spatial profiles of the magnitude of the total EF and of its components 
(Fig. 6). Grid montages are more similar to 2-electrode montages with inter- 
electrode distance of 5 mm, compared to the montages with 65 mm. The latter also 
originate EF magnitude larger than 0.50 Emax in a wider region, by approximately 
two spinal segments than the grid or the 2-electrode at 5 mm. 
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Fig.6 Magnitude of the EF and components in spinal-WM averaged over 1-mm thick slices along 
the spinal cord length and normalized to EF maximum (Emax). The corresponding electrode mon- 
tage and orientation and the legend for each element are indicated in each plot at the top left and 
right, respectively. The E; component is superimposed on the total EF magnitude in the first top row 


The direction of the EF is consistent with the relative position of anodes and 
cathodes in each case. TO montages originate an EF with a E, component that 
almost matches the total magnitude (Fig. 6, top row). The diagonal and longitudinal 
alignments (D45 and L90) increase E,,,. and E,, contributions to the total EF (Fig. 6, 
middle and bottom rows). E,, contribution is similar to the other components in D45, 
decreasing more in 2-D at 65 mm; this contribution disappears almost entirely at 
L90 montages. 

The Eng component peaks in the anode-cathode transition region; E,, presents 
two peaks of opposite sign, one below each anode and cathode regions. The larger 
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distance between electrodes in 2-D and 2-L at 65 mm increases E,, below the elec- 
trodes, resulting in a total EF spatial profile with two peaks, below anode and cath- 
ode regions. The anode-cathode relative placement and distance have a strong 
influence on the direction of the EF, allowing to establish which component can be 
more significant to the total EF magnitude. This could be determinant to target spe- 
cific spinal neurons according to their orientation inside the SC. 


3.4 Semi-homogeneous Models: How Tissues’ Different 
Conductivities Influence the EF 


The EF spatial profiles in the semi-homogeneous models can reveal the influence of 
the electrical conductivity of each tissue on EF components and total magnitude 
along the SC. The EF spatial profile is very similar to the profile for the Ej; compo- 
nent in TO montages in the spinal WM (Fig. 7). The same type of profiles occurs in 
the spinal GM (not shown). This profile presents peak-like features in the spinal 
segments below and near the Grid or the 2-electrode. Also, the total magnitude and 
E, profiles mimic the semi-homogeneous vertebrae model profile, thus the conduc- 
tivity of vertebrae seems to be the main factor determining the existence of the 
peak-like features of the EF spatial variation in the spinal WM (Fig. 7, light grey 
dashed line in each plot). The E,, component contributes to 97-99% of Emax in all TO 
models, which highlights the importance of the anode-cathode placement in the 
total EF spatial distribution (Fig. 10). 

D45 montages show similar profiles for the total EF and components (Fig. 8). It 
presents the same type of global profiles as in the homogeneous model, where Ej, 
represents 69-75% of Emax in the homogeneous models. This percentage is around 
95-99% for the heterogeneous (isotropic) model, the larger value corresponding to 
2-D at d = 65 mm (Fig. 10). This increase in Ej, contribution should arise due to 
the influence of the CSF conductivity, since the total EF has a profile more similar 
with semi homogeneous CSF model (Fig. 8, dot line in dark grey). There are also 
some peak-like features in the spatial distribution at the same spinal level as in the 
semi-homogeneous vertebrae model, which supports the influence of the vertebrae 
conductivity in local hotspots. The E, component is larger in the semi-homogeneous 
models, presenting the same peaks below the anode and cathode regions as in the 
homogeneous model, but almost disappears in the heterogeneous isotropic model. 
This indicates that the combined effect of the different conductivities cancels out the 
contribution of the electrode placement to the E,, component. 

L90 montages have similar distributions to D45 montages, with the Ej, compo- 
nent contributing to 99% of Emax value, and with an almost negligible E; component 
(Figs. 9 and 10). The E,, is present in homogeneous and semi-homogeneous mod- 
els, but the overall effect of different conductivities results in a less meaningful 
contribution in the heterogeneous model, just as seen in the D45 montages. 
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Fig.7 EF profiles of the total and components magnitude normalized to Emax in the spinal-WM for 
TO montages. The first two rows present the spatial profiles for the total EF, Eiong, Eva and E; for 
Grid 0? montage. The middle and bottom rows present the spatial profiles for the total EF and E; 


for Transversal 2-pads at d = 5 mm and d = 65 mm, respectively 
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z (mm) 


Fig. 11 Distribution of vertebral volume, disk volume and EF magnitude in Grid 0? montage, 
averaged over 1-mm thick slices along the SC length and normalized to maximum values in each 
case (horizontal bottom and upper axes correspond to EF and vertebral/disk volume, respectively). 
To the left, the EF magnitude volume distribution is represented with the same length for compari- 
son and uses the same colour scale as in Fig. 4. The red lines show the spatial correspondence 
between local hotspots, EF peaks, disk volume maxima and vertebral volume minima 


What determines the “peak-like” profile seen in the E, component? Figure 11 
compares the normalized EF magnitude profile with the distribution of volume of 
vertebrae and disks for Grid 0? montage, where the E,, component largely influences 
the total EF magnitude. The positions of vertebral spaces and intervertebral disks 
(volume minima and volume maxima, respectively) are coincident with the “peak- 
like" features in the EF spatial distribution. Considering that the profile of semi- 
homogeneous vertebrae is the only one that reflects these “peak-like” features, 
vertebral conductivity must have the largest influence in the formation of local 
hotspots. 
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4 Discussion 


This study addresses the interplay of different electrode placements and electrical 
conductivity modelling paradigms in setting the main characteristics of the EF 
induced by tsDCS over the spinal cord. Our main finding is that the relation between 
the two factors is complex: anode-cathode placement determines the EF orientation, 
however electrical properties of tissues can change this orientation and originate 
local hotspots. Conversely, local hotspots that could be originated by tissue conduc- 
tivity heterogeneities can disappear when the electrodes are oriented along the SC, 
by increasing the contribution of the highly conductive CSF. This variability in the 
interplay between electrodes and electrical properties of tissues in the current path 
is determinant to optimize EF delivery using tsDCS for a specific spinal target. 


4.1 Anode-Cathode Placement Interplays with Tissues’ 
Conductivities to Define the EF Direction 


The homogeneous model is a useful tool to isolate the effect of electrode placement 
since it considers the human model as a completely uniform volume conductor. The 
main observations using this model are: 


* Different electrode geometries and number are not determinant in the EF direc- 
tion if the anode-cathode spatial relation is preserved — a grid of 16 small circular 
electrodes or a two-electrode configuration resulted in similar EF spatial profiles, 
if presenting the same anode-cathode relative position (Fig. 6); 

* Larger distances between electrodes originate higher EF magnitudes in wider 
regions (Fig. 6, Table 2). 


When turning to the isotropic heterogeneous model, the EF shows peak-like fea- 
tures and does not always preserve the spatial profiles of the EF and of its compo- 
nents (Figs. 4, 7, 8, and 9). For instance, E,, profile in D45 and L90 homogeneous 
models has two main peaks and with opposite orientations below anodes and cath- 
odes, that almost disappear in the heterogeneous model. This indicates that electri- 
cal conductivities considered in the heterogeneous (isotropic) model interfere with 
the EF orientation previously established by anode-cathode relative position. This 
raises the question: which tissues have the most impact on the EF spatial distribu- 
tion in the spinal cord? Realistic numerical models of DC cortical stimulation refer 
CSF and skull as the tissues that can induce large variability in EF magnitude [13, 
14]. Modelling studies on tsDCS revealed a negative relation between CSF volume 
and EF magnitude over the thoracic region when two electrodes are placed over T10 
spinous process and right arm [7]. In our previous modelling study on thoraco- 
lumbar tsDCS, locations of vertebrae bony edges and CSF narrowing were associ- 
ated with local hotspots, regardless of electrode placement [9]. The 
semi-homogeneous models were considered above to isolate the contribution of 
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each tissue’s conductivity. In all cases, vertebrae and CSF are the two tissues that 
have the most impact in local features and contribute to the global spatial profile. 
The EF profiles in vertebrae and CSF semi-homogeneous models closely resemble 
the EF in the heterogeneous model for TO and D45/L90 montages, respectively. 
Vertebrae and CSF conductivity values differ by almost one order of magnitude 
relative to those of other tissues, which may explain a larger effect in current spread- 
ing to the surrounding tissues, thus originating the profiles observed: vertebrae has 
alow conductivity, which favors the passage of current through inter-vertebral space 
to the spinal cord, favoring the E4 component; CSF’s higher conductivity contrib- 
utes to an increase in the longitudinal component of the EF and also contributes to 
current focusing effects in narrower regions originating, for instance, the sharp peak 
observed at T12 level in D45/L90 montages. The effect of CSF thinning was also 
observed in transcranial magnetic stimulation models, where the induced current 
flows tangential to the skull, leading to a local increase of the EF magnitude at the 
gyral crowns due to local CSF thinning [22]. 

The effect of distance was already addressed in previous tsDCS modelling works. 
Anode-cathode relative positions and tissues conductivities should be carefully con- 
sidered when choosing a tsDCS montage. The interplay between these two factors 
will determine the EF orientation. The modulation of spinal neurons, just as in corti- 
cal neurons, will occur if these are aligned with the induced EF. Cellular models of 
spinal motor neurons identified axon terminals as the dominant cellular target of 
tsDCS and misplacements of 5 cm in electrodes’ positions lead to a change in EF 
direction and a consequent reversal of polarization at target, thus the influence of 
electrode distance in position should be carefully considered [8]. However, DCS 
does not change excitability of peripheral axons, which suggests that synapses are 
the main targets for excitability modulation induced by DCS [23, 24]. Spinal dys- 
function can have many causes but are most frequently associated with aberrant 
triggering of spinal reflexes and muscle tone, due to vertebral lesions, tumors or 
neuron degeneration. Many different neurodegenerative, inflammatory or traumatic 
disorders can affect spinal cord networks, with variable modulation by the same 
tsDCS montage. This leads to some considerations when translating modelling find- 
ings to the clinical context: 


1. The selection of the topography and polarity of electrodes in tsDCS is strongly 
dependent on the position and orientation of the targeted spinal neurons position 
and their nearby anatomical frame — it is important to consider the location of the 
surrounding vertebrae and nearby CSF narrowing due to protruded or herni- 
ated disks; 

2. Electrode grid should position the target below the anode-cathode transition 
region; two-electrode montages should comprise the target between the posi- 
tions of the two electrodes; 

3. Longitudinal/diagonal placements will be more favorable to modulate 
longitudinally-oriented targets; transversal placements will preferentially modu- 
late neural targets with a mediolateral orientation; 

4. In longitudinal/diagonal placements, increasing the distance between the elec- 
trodes increases the EF magnitude at the cost of focality. 
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4.2 Multi-electrode Montages Can Be Relevant for Differential 
Targeting in tsDCS 


Safety concerns are a major issue when considering transcutaneous application of 
direct currents. Most applications of tsDCS have been more conservative, with deliv- 
ery of current intensities below 3 mA. However, intensities up to 4 mA are being 
considered to increase neuromodulation outcomes in brain stimulation, with promis- 
ing results regarding patients' tolerability [21]. The use of a multi-electrode grid may 
raise some issues because current density will be larger than 2-electrode montages at 
the skin located below the anode-cathode transition region. However, our calculations 
predict maximum values around one order of magnitude below the threshold for tissue 
damage [19]. Also, multi-electrode grid montages can reproduce the same type of EF 
spatial distribution than the traditional two-electrode montages, as demonstrated in 
this study, however with lower EF values. The maximum EF obtained in the grid rota- 
tions considered above are larger than 0.15 V/m in the heterogeneous (isotropic) 
model (Emax (grid 0?) = 0.166 V/m; Emas (grid 45°) = 0.320 V/m; Ema (grid 
90°) = 0.438 V/m; Fig. 4), an assumed threshold value for effective neuromodulation 
(Sect. 2.3). What would be the advantage of multi-electrode grid montages, if the EF 
orientations reproduced with the rotations can be obtained with 2-electrode montages 
and with higher EF magnitude? In the clinical setting, a multi-electrode grid system 
may be an advantage if spinal targets with different orientations are considered, 
because the grid allows to obtain different anode-cathode alignments without chang- 
ing the placement of the grid, for example, to have a stimulation paradigm alternating 
mediolateral and longitudinally oriented EFs. A grid setting using only 8 of the 16 
electrodes can reproduce an EF with a spatial distribution similar to the D45 group, by 
considering the 4 top right electrodes as anodes and the 4 bottom left electrodes as 
cathodes and leave the other electrodes turned off (Fig. 12). Thus, the multi-electrode 


EF magnitude EF magnitude 


— EF longitudinal 


1650 | — EF longitudinal 


€t 

Li 

cs 

G ' 
cs 1600 4 9 

EF right-left Co , 

c 

€x 

n 

T 


— EF wentral-donal = EF ventral-dorsal 


EF right-left 


71 
E 
eo. 
os 


ee 
150 | "9 


— spinal segments 


— spinal segments 


i Emax = 0.460 V/m 


z (mm) 


EF/EF,... in spinal-WM 


Fig. 12 Distribution of total EF magnitude and EF components in the diagonal grid (left) and grid 
rotation 45? (right) montages, averaged over 1-mm thick slices along the SC length and normalized 
to maximum values in each case. The EF magnitude volume distribution over T6-Co segments in 
WM is represented at the top left corner of each plot, with the value of Emax 


120 S. R. Fernandes et al. 


grid can be used for differential targeting, i.e. to stimulate different spinal networks 
according to their orientation at specific periods within the same session, which can be 
of great relevance for clinical applications of tsDCS for multi-factorial spinal 
dysfunctions. 


4.3 Limitations in Modelling tsDCS: What Lies Ahead 


The present study only considers isotropic electrical conductivities. In a previous 
study, we included artificial anisotropy tensors of WM and muscle, considering the 
typical direction of fibers. Muscle anisotropy only contributed to change slightly the 
EF magnitude, however the anisotropy of WM increased the EF sensitivity to ana- 
tomical characteristics and spinal curvature, which can have an impact on the pre- 
dicted neuromodulation on certain spinal regions. Thus, we recommend that the use 
of models to optimize tsDCS protocols should include, whenever possible, anisotro- 
pic considerations for the WM and GM. 

The model presented here has an artificial design of the GM and does not repre- 
sent the spinal rootlets with dorsal ganglia, where sensory neurons somas are 
located. Realistic representation of these components will be a difficult step, since 
it will require MRIs of very high resolution, with a segmentation procedure which 
will take many semi-automatic or even manual time-consuming tasks. Even so, this 
will be a necessary update of human realistic models to address the true neuromodu- 
latory potential of tsDCS. 

Connecting macroscale representations of the SC with neuronal and circuit mod- 
els is also essential to determine what are the effects at neuronal and synaptic level. 
The long-term effects of tsDCS may be similar to the processes of long-term poten- 
tiation (LTP) observed in cortical stimulation, related with neuroplasticity [18, 25]. 
Neuronal and network models can unveil how tsDCS may be applied to repair spi- 
nal network communication and recovery of function at long-term, thus providing a 
solid support of tsDCS as an effective therapeutic resource for sensorimotor 
rehabilitation. 


5 Conclusion 


Non-invasive stimulation of the spinal cord in the form of tsDCS has a promising 
therapeutic potential to address sensorimotor dysfunctions of the spinal cord. 
Computational numerical studies can provide useful information to optimize the 
delivery of currents, indicating how to combine knowledge of the electrical proper- 
ties of tissues and relative placement of electrodes to increase stimulation selectivity 
aiming at a specific spinal target. Even if not very different from traditional 
2-electrode montages, multi-electrode grid stimulation paradigms can provide dif- 
ferential stimulation, by changing EF directions over different periods without 


Interplay Between Electrical Conductivity of Tissues and Position of Electrodes... 121 


changing the position of the electrode grid. This could be helpful to modulate dif- 
ferent spinal circuits in a single session, when these vary in orientation inside the 
spinal cord. Furthermore, the grid system may be further optimized by adjusting 
electrode size and inter-electrode distance, since these characteristics can change 
EF magnitudes, as predicted in 2-electrode montages. 
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Part V 

High Frequency Electromagnetic Modeling 
and Experiment: MRI Safety with Active 
and Passive Implants 


RF-induced Heating Near Active A 
Implanted Medical Devices in MRI: gese 
Impact of Tissue Simulating Medium 


James E. Brown, Paul J. Stadnik, Jeffrey A. Von Arx, and Dirk Muessig 


1 Introduction 


Generally, patients with active implanted medical devices (AIMD) are denied 
access to MRI due to the potential for hazardous interactions. It has been estimated 
that 17% of pacemaker patients will need an MRI within 12 months of device 
implantation [1]. MR imaging is a highly desirable imaging modality due to its non- 
ionizing radiation and image quality for soft tissue imaging [2]. Therefore, in recent 
years device manufacturers, academics, and regulatory agencies have developed 
standardized test methods for improving the access to this important diagnostic tool 
for AIMD patients [3, 4]. 

Several potentially hazardous interactions of implanted medical devices with the 
MR system have been identified [5, 6]. Computational human models (CHMs) are 
used in the evaluation of three hazards: RF-induced heating, RF-induced unintended 
stimulation, and RF-induced malfunction. 


2 Utilizing Computational Human Models 
for the Assessment RF-induced Heating 


Computational human models have been used for a range of electromagnetic appli- 
cations [7]. This chapter discusses the specific application of CHMs to MRI 
RF-induced heating near AIMDs. While some studies have used measurement 
methods in homogenous phantoms [8], this work focuses on the methods outlined 
in [3, 4], which are based on the established transfer function method [9, 10]. The 
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benefits of using CHMs in this work have been outlined in [11, 12], which provide 
a framework for evaluating millions of scenarios and can be contrasted with the test- 
ing of pacemaker resilience to Electronic Article Surveillance (EAS) systems [13], 
which involves physically testing the device in a more limited number of prescribed 
tests. The general process is illustrated in Fig. 1. 

The standardized test methods of [3] outline a 4-tier system for estimating in vivo 
RF-induced heating utilizing CHMs. Lower-numbered tiers are more conservative 
than higher-numbered tiers, that is Tier 1 is the most conservative and Tier 4 is con- 
sidered the most accurate. Tier 1 is impractical and will be removed in the upcoming 
standard [14]. 

Tiers 2 and 3 both involve simulating a set of CHMs and a set of MR birdcage 
coils to build a library of expected field distributions within the human body. 
Examples of CHMs include the Virtual Family [15] and the Visible Human Project 
[16]. By varying parameters such as landmark position (alignment of the center of 
the MR coil with the CHM), body position, coil geometry, and using a variety of 
CHMs (different heights, weights, BMI, ages, efc.), millions of scenarios can be 
investigated. In a Tier 2 analysis, the risk to the patient is assessed using the maxi- 
mum fields over the expected implant volume, after taking a 10 g average. In a labo- 
ratory setting, the device is then exposed to the appropriate field level and the 
temperature rise is measured. For elongated devices such as leaded AIMDs, as the 
maximum field values are taken over a large volume and then applied uniformly 
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Model E-field distribution in multiple Generate AIMD System Model 
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models for each AIMD system configuration 
(device / lead combination) 


human body model variants, RF coils 
and positions. 


Calculate and Analyze Deposited Power at Lead Tip 
Calculate deposited power at tip for multiple pathways in each human body model variant 


for each AIMD system configuration. 


Determine Animal Study Target 
Determine target deposited power from the distribution 


for each lead in the system. Scale by Uncertainty. 


Animal Study 
Confirm safety for each lead (based on lead power deposition at or 


above animal study target). 


Fig. 1 Using CHMs to assess risk to the patient due to RF-induced tissue heating 
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over the entire device, this would result in impractically high test conditions and 
severely over-predict RF-induced heating. A Tier 3 analysis uses the same library of 
3-D field distributions, but extracts the incident fields over one or more lead path- 
ways and then combines this with the separately-derived lead model. Then, a prob- 
ability distribution function of the estimated in vivo RF-induced heating can be 
carried forward to a safety assessment using an animal model or alternate method of 
clinical assessment. 

The lead models used in Tier 3 analyses are often piecewise excitation models, 
where the lead response (e.g., RF-induced deposited power or temperature rise near 
the lead tip) is quantified for the excitation of a section of the lead (e.g., 1 cm). The 
so-called reciprocity method [10] involves a single excitation and a swept measure- 
ment of the response. Whether numerically- or experimentally-derived [17-21], 
transfer function models are developed within a homogenous medium, commonly 
referred to as a tissue-simulating medium (TSM). This work focuses on the applica- 
tion of the homogenous TSM used in the lead model to the inhomogeneous CHM 
used in the generation of the incident fields, and the ultimate impact to the in vivo 
estimation of RF-induced heating. 


3 Tissue-Simulating Media (TSM) 


Equivalent medium theorem [22] suggests that there exists an optimal homogenous 
TSM for agreement between the lead model (transfer function) and predicted in 
vivo RF-induced heating (in inhomogeneous CHMs). Indeed, recommended proce- 
dures for identifying such a medium are becoming standard practice [13] for AIMD 
manufacturers. The following subsections will discuss what is known about the 
impact of TSM on the AIMD model. 


3.1 Effect of TSM on Computation of RF-induced Heating 


The analytical problem of a medical device lead exposed to RF fields from an MR 
scanner is a particular instance of a wire scatterer in an incident field. Classically, 
these problems are investigated via the Method of Moments (MoM) [23], though 
MoM is known to be disadvantageous for analysing problems with multiple mate- 
rial boundaries as are present in the human body. Further, the problem of bare and 
insulated wire (as antennas or scatterers) in an extended conductive medium has 
been extensively studied in the literature [24—28], especially for submarine com- 
munication and geophysical exploration. 

This analytical framework has been extended to the particular problem of 
implanted medical device leads [29], whereby the particular impact of the TSM to 
the predicted in vivo RF-induced heating can be quantified. Using a safety index 
which is proportional to the elevated SAR near the lead tip, general trends are seen 
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for the impact of the TSM on resonant characteristics of the lead geometry. The 
resonant effect has also been extensively defined in the literature [30-32]. One par- 
ticularly important consequence to the assessment of RF-induced heating near the 
lead tip is the necessity of elevated Electric field at the lead tip due to the discontinu- 
ity in the current distribution, which can be shown using a basic MoM formulation 
[33]. While the preceding analytical techniques focus on quantifying RF-induced 
heating in terms of electromagnetic quantities, additional considerations of the 
TSM for measuring temperature rise are also required [34, 35]. 

While these studies provide an understanding of the effect of changing the sur- 
rounding medium, or TSM, on the lead response, more contemporary research has 
focused on the impact of the choice of TSM on the overall accuracy of the in vivo 
prediction of RF-induced heating. This important work is discussed in the next two 
subsections, which are divided into studies of simplified structures such as a simple 
wire and more complex geometries representing actual geometries of practical 
AIMD leads. 


3.2 Numerical Studies of Simplified Structures 


A basic model for an AIMD lead is an insulated wire with a bare section at the distal 
end. This simplified structure is useful for its ease in construction and for simplify- 
ing analyses, and thus is used for standardization of test methods [3]. Indeed, such 
structures are often used for numerical simulations of transfer functions, examples 
include [18, 19]. 

A transmission-line model of a simplified structure was shown in [36] to give 
good agreement with measurement for a variety of parameters of the wire (insula- 
tion thickness, conductor diameter, insulator dielectric constant, lead length, etc.) 
and for different TSMs, for RF-induced voltage at the proximal end. While this 
work did not focus on RF-induced heating, it is still important for any discussion of 
the impact of TSM properties to highlight the influence of the tip impedance, Z,;,, on 
the ultimate result, a relationship which has been explored in further depth [37]. 
Still, neither of these works investigated which TSM would then give the most accu- 
rate result when compared to a full evaluation with an inhomogeneous CHM. This 
question for simplified structures is answered in [38], which concludes that the 
appropriate choice in homogeneous TSM is the TSM with the conductivity which 
matches that at the lead tip. This conclusion is dependent on the thickness of the 
insulation surrounding the conductor, as [39] shows that for a given TSM the SAR 
will plateau for increased insulation thickness. For thin insulation leads, the TSM at 
the lead tip is no longer dominant and the choice in TSM should align with the aver- 
age along the lead body [40]. While the focus of this work has been RF-induced 
heating, most of these principles are extensible to other RF-induced hazards, such 
as RF-induced unintended stimulation [41]. 
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3.3 Numerical Studies of Realistic Device Geometries 


Medical device leads are of course more complex than the simplified structures of 
the previous subsection. For example, pacemaker leads are helically-wound, with 
multiple conductors present which may be coaxial or coradial [42]. The current ver- 
sion of the ISO technical specification [3] for MRI safety of AIMDs includes a set 
of suggested TSMs for manufacturers to consider. For electrically short, non-leaded 
devices the analysis may be more straightforward. An example of this analysis is 
shown in [43], which follows the test methods defined in [4]. In this type of analysis, 
typically two TSMs are used, and the worst-case of these two scenarios is carried 
forward. In some situations, this may be overly conservative. 

Importantly, the set of geometries in [44] included up to four helically wound 
conductors and includes a parametric study of the helix design, to explore the rela- 
tionship of the tip impedance to the predicted RF-induced heating. The importance 
of the lead pathway to the optimal choice of TSM is discussed, and a novel method 
for constructing an inhomogeneous computational phantom for transfer function 
generation is proposed. 

For a percutaneous SCS lead, the choice of the high conductivity medium (HCM) 
from [3] was shown to be more accurate, while the low conductivity medium (LCM) 
would lead to an over-prediction of RF-induced heating by 72-74% [45]. This con- 
trasts with another study of a simplified pacemaker lead [46] which shows that 
LCM is more accurate while the HCM overpredicted RF-induced heating by 80%. 
However, exact insulation thicknesses for the lead models used in these two studies 
were not given. Both studies illustrate the importance of TSM selection on the pre- 
dicted risk to the patient due to RF-induced heating. While safety assessments 
should always be conservative, care must be taken to avoid situations where the 
prediction is so high as to prevent patient access to this important diagnostic tool. 


4 Discussion 


In this chapter, we have investigated the impact of the CHM on the estimation of 
MRI RF-induced heating near the tips of AIMD leads through the lens of TSM 
selection during transfer function development. While general guidance for TSM 
selection is still being investigated, early evidence exists to recommend that device 
manufacturers use a TSM similar to that found at the distal end for leads with thick 
insulation, while for thin insulation, a homogenous TSM should approximate the 
average properties of tissues in contact with the lead body. The data is supported by 
experimental, theoretical, and computational analysis of simplified structures as 
well as practical lead geometries. 
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41 Future Work 


The techniques shown in this work are currently being extended to provide general 
guidelines for the selection of TSM and its applicability to the CHM. AIMD manu- 
facturers desire to have a more defined transition from "thick insulation" to “thin 
insulation", a definition which may vary due to the particular lead construction 
(straight versus helical wires, number of conductors, diameter of the conductors 
versus outer diameter of the lead, efc.). In addition, these types of sensitivity analy- 
ses are enabled mostly by numerical simulation, as to parameterize physical quanti- 
ties could make prototyping and testing time-consuming and expensive. 

The further refinement of numerically derived transfer functions can also enable 
more solutions for mitigating MRI RF-induced heating. Current practice is still to 
build prototypes and experimentally confirm the reduction in expected RF-induced 
heatingRF-induced heating [47, 48]. The investigation of solutions can only be 
improved by a more thorough understanding of the influence of the computational 
model on the ultimate results of the safety assessment. 


References 


1. R. Kalin, M.S. Stanton, Current clinical issues for MRI scanning of pacemaker and defibrilla- 
tor patients. Pacing Clin. Electrophysiol. 28(4), 326-328 (2005) 

2. M.S. Brown, R.C. Semelka, Concept of magnetic resonance, in MRI: Basic Principles and 
Applications, 3rd edn., (John Wiley & Sons, Hoboken, 2004), pp. 11-20 

3. ISO/TS 10974:2018 (E), "Assessment of the safety of magnetic resonance imaging for patients 
with an active implantable medical device", 2018 

4. ANSI/AAMI PC76:2021, Requirements and Test Protocols for Safety of Patients with 
Pacemakers and ICDs Exposed to MRI, 2021 

5. L.P. Panych, B. Madore, The physices of MRI safety. J. Magn. Reson. Imaging 47, 28—43 (2018) 

6. J. Kabil et al., A review of numerical simulation and analytical modeling for medical devices 
safety in MRI. Yearb. Med. Inform., 152-158 (2016) 

7. S.N. Makarov et al., Virtual human models for electromagnetic studies and their applications. 
IEEE Rev. Biomed. Eng. 10, 95-121 (2017) 

8. P. Nordbeck et al, Spatial distribution of RF-induced E-fields and implant heating in 
MRI. Magn. Reson. Med. 60, 312—319 (2008) 

9. S.-M. Park, R. Kamondetdacha, J.A. Nyenhuis, Calculation of MRI-induced heating of an 
implanted medical Lead wire with an electric field transfer function. J. Magn. Reson. Imag. 
26, 1278-1285 (2007) 

10. S. Feng et al., A technique to evaluate MRI-induced electric fields at the ends of practical 
implanted Lead. IEEE Trans. Microw. Theory Techn. 63(1), 305-313 (2015) 

11. E. Brown, et al., MR Conditional Safety Assessment of Implanted Medical Devices: 
Advantages of Computational Human Phantoms, Proc. 36th Annu. Int. Conf. IEEE EMBC, 
Orlando, FL, pp. 6465-6468, 2016 

12. B.L. Wilkoff et al., Safe magnetic resonance imaging scanning of patients with cardiac rhythm 
devices: a role for computer modeling. Heart Rhythm. 10(12), 1815-1821 (2013) 

13. R. Herkert, E3 Test Protocol for Medical Devices to Security and Logistical Systems, Ver. 6.0 
(Medical Device Test Center, Georgia Tech Research Institute, Atlanta, 2013) 


RF-induced Heating Near Active Implanted Medical Devices in MRI: Impact of Tissue... 131 


14. 


15. 


16. 


17; 


18. 


19. 


20. 


21. 


22. 


23. 
24. 


25. 


26. 


27. 


28. 


29. 


30. 


31. 


32; 


33. 


34. 


35. 


36. 


ISO 10974 (Draft), Assessment of the safety of magnetic resonance imaging for patients with 
an active implantable medical device, To be published 

A. Christ et al., The virtual family—development of surface-based anatomical models of two 
adults and two children for dosimetric simulations. Phys. Med. Bio. 55, N23-N38 (2010) 
G.M. Noetscher, et al., Computational Human Model VHP-Female Derived from Datasets 
of the National Library of Medicine, Proc. 38th Annu. Int. Conf. IEEE EMBC, Orlando, FL, 
2016, pp. 3350-3353 

A. Yao et al., Efficient and reliable assessment of the maximum local tissue temperature 
increase at the electrodes of medical implants under MRI exposure. Bioelectromagnetics 
40(6), 422—433 (2019) 

M. Kozlov, W. Kainz, Comparison of lead electromagnetic model and 3D EM results for helix 
and straight leads, Proc. 19th Int. Conf. Electromagn. Adv. Appl., pp. 649-652, 2017 

M. Kozlov, W. Kainz, Lead electromagnetic model to evaluate RF-induced heating of a coax 
lead: a numerical case study at 128 MHz. IEEE J. Electromagn. RF Microw. Med. Biol. 2(4), 
286—293 (2018) 

E. Zastrow, M. Capstick, N. Kuster, Experimental system for RF-heating characterization of 
medical implants during MRI’, Proc. 24th Annu. Meeting ISMRM, Singapore, 2016 

E. Zastrow, A. Yao, N. Kuster, Practical considerations in experimental evaluations of 
RF-induced heating of leaded implants, 32nd URSI GASS, Montreal, Canada, 2017 

Y. Wang et al., On the development of equivalent medium for active implantable device radio- 
frequency safety assessment. Magn. Reson. Med. 82, 1164—1176 (2019) 

R.F. Harrington, Field Computation by Moment Methods (IEEE Press, New York, 1993) 

P.E. Atlamazoglou, N.K. Uzunoglu, A Galerkin moment method for the analysis of an insulated 
antenna in dissipative dielectric medium. IEEE Trans. Microw. Theory 46, 988—996 (1998) 
R.W.P. King, G.S. Smith, Antennas in Matter: Fundamentals, Theory, and Applications (The 
MIT Press, Cambridge, MA, 1981) 

R.W.P. King, B.S. Trembly, J.S. Strohbehn, The electromagnetic field of an insulated antenna 
in a conducting or dielectric medium. IEEE Trans. Microw. Theory Techn. MTT-31(7), 
574-583 (1983) 

R.W.P. King, Antennas in material media near boundaries with application to communica- 
tion and geophysical exploration, part I: the bare metal dipole. IEEE Trans. Ant. Propagat. 
AP-34(4), 483-489 (1986) 

R.W.P. King, Antennas in material media near boundaries with application to communica- 
tion and geophysical exploration, part II: the terminated insulated antenna. IEEE Trans. Ant. 
Propagat. AP-34(4), 490—496 (1986) 

C.J. Yeung, R.C. Susil, E. Atalar, RF safety of wires in interventional MRI: using a safety 
index. Magn. Reson. Med. 47, 187-193 (2002) 

J.E. Brown, C.S. Lee, Radiofrequency resonance heating near medical devices in magnetic 
resonance imaging. Microwave Opt. Technol. Lett. 55(2), 299—302 (2013) 

S.O. McCabe, J.B. Scott, Cause and amelioration of MRI-induced heating through medical 
implant lead wires, 2/st Elect New Zealand Conference, Hamilton, New Zealand, Nov 2014 
S.O. McCabe, J.B. Scott, Technique to Assess the Compatibility of Medical Implants to the RF 
Field in MRI, Asia-Pacific Microwave Conference 2015 6—9 Dec 2015 

J.E. Brown, Radiofrequency heating near medical devices in magnetic resonance imaging, 
Ph.D. dissertation, Bobby B. Lyle School of Engineering, Southern Methodist University, 
Dallas, TX, 2012 

S.M. Park et al., Gelled vs. non-gelled phantom material for measurement ofMRI-induced 
temperature increases with bioimplants. IEEE Trans. Magn. 39(5), 3367—3369 (2003) 

C.D. Smith, J.A. Nyenhuis, K.S. Foster, A comparison of phantom materials used in evaluation 
of radiofrequency heating of implanted medical devices during MRI, Proc. 23rd Annu. Int. 
Conf. IEEE EMBC, Istanbul, Turkey, pp. 2311-2314, 2001 

J. Liu et al., A transmission line model for the evaluation of MRI RF induced fields on active 
implantable medical devices. IEEE Trans. Microw. Theory Techn. 66(9), 4271-4281 (2018) 


132 J. E. Brown et al. 


37. J. Liu, et al., On the relationship between impedances of active implantable medical devices 
and device safety under MRI RF emission, IEEE Trans. EMC, 2019 (Early Access) 

38. K.N. Kurpad, et al., MRI RF safety of Active Implantable Medical Devices (AIMDs): numeri- 
cal study of the effect of conductivity of tissue simulating media on device model accuracy, 
Proc. 26th Annu. Meeting Int. Soc. of Magn. Reson. Med., Paris, France, pp. 4075, 2018 

39. P.A. Bottomley et al., Designing passive MRI-safe implantable conducting leads with elec- 
trodes. Med. Phys. 37(7), 3828—3843 (2010) 

40. J.E. Brown, et al., MRI safety of active implantable medical devices: numerical study of the 
effect of lead insulation thickness on the RF-induced tissue heating at the lead electrode, 43rd 
Annu. Int. Conf. IEEE EMBC, 2021 

41. J.E. Brown et al., RF-induced unintended stimulation for implantable medical devices in 
MRI, in Brain and Human Body Modeling 2020: Computational Human Models Presented 
at EMBC 2019 and the BRAIN Initiative® 2019 Meeting, ed. by S. Makarov et al., (Springer 
Nature, Cham, Switzerland), pp. 283-292 

42. C. Tang etal., Initial experience with a co-radial bipolar pacing lead. Pacing Clin. Electrophysiol 
20(7), 1800-1807 (1997) 

43. J.E. Brown et al., Calculation of MRI RF-induced voltages for implanted medical devices using 
computational human models, in Brain and Human Body Modeling: Computational Human 
Modeling at EMBC 2018, ed. by S. Makarov et al., (Springer Nature, Cham, Switzerland), 
pp. 283-294 

44. J. Liu et al., Investigations on tissue-simulating medium for MRI RF safety assessment for 
patients with active implantable medical devices. IEEE Trans. EMC 61(4), 1091-1097 (2019) 

45. X. Min, S. Sison, Transfer functions of a spinal cord stimulation systems in mixed media and 
homogeneous media for estimation of RF heating during MRI scans, Proc. 40th Annu. Int. 
Conf. IEEE EMBC, Honolulu, HI, 2018, pp. 2048-2051 

46. X. Min, S. Sison, Impact of mixed media on transfer functions with a pacemaker system for 
estimation of RF heating during MRI scans. Comput. Cardiol. 44, 1—4 (2017) 

47. P. Nordbeck et al., Reducing RF-related heating of cardiac pacemaker leads in MRI: imple- 
mentation and experimental verification of practical design changes. Magn. Reson. Med. 68, 
1963-1972 (2012) 

48. P. Serano et al., Novel Brain stimulation technology provides compatibility with MRI. Sci. 
Rep. 5, 9805 (2015) 


Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0 
International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, 
adaptation, distribution and reproduction in any medium or format, as long as you give appropriate 
credit to the original author(s) and the source, provide a link to the Creative Commons license and 
indicate if changes were made. 

The images or other third party material in this chapter are included in the chapter's Creative 
Commons license, unless indicated otherwise in a credit line to the material. If material is not 
included in the chapter's Creative Commons license and your intended use is not permitted by 
statutory regulation or exceeds the permitted use, you will need to obtain permission directly from 
the copyright holder. 


Computational Tool Comprising Visible (R) 
Human Project® Based Anatomical M 
Female CAD Model and Ansys HFSS/ 
Mechanical& FEM Software 

for Temperature Rise Prediction Near 

an Orthopedic Femoral Nail Implant 
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1 Introduction 


This tool partially described in Refs. [1—4] is a non-clinical assessment model used 
to predict the RF power deposition induced temperature rise near orthopedic 
implants as a function of implant geometry, location, material, and scan time. It is 
using the complete in silico multiphysics MRI environment and the anatomical 
human CAD model. The validated context of use is limited to a mid-aged or elderly 
female subject of 50—70 years old with a higher obesity (BMI or body mass index 
of 30-36) scanned in a 1.5 T full-body circularly polarized cylindrical MRI bird- 
cage coil with two ports at 64 MHz. It is also limited to one yet critical implant 
geometry: a long femoral nail (a nearly straight long metal rod) which is subject to 
most excessive heating during long scan times at 1.5 T. 

This tool can augment the widely used ASTM F2182 standard that measures RF 
implant heating in a homogeneous gel phantom by providing extra safety margins 
caused by the influence of the realistic heterogeneous human body. It can help to 
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identify the appropriate worst-case implant size, configuration, and orientation as a 
function of the scan protocol and required scan time (the output RF power). The tool 
does provide the estimates of SAR values in the regions around the implant. The 
tool does not provide heating or SAR estimates for gradient-coil induced power 
deposition for orthopedic implants. 

A number of modeling and physical processes were tested to validate this tool. 
They include: (i) topological validation of the entire female CAD model (required 
approximately 4000 man hours for its construction); (ii) anatomical validation of 
the constructed human model by anatomical experts; (iii) validation of the human 
model and FEM software compatibility; (iv) SAR deposition validation in the 1.5 T 
full-body MRI birdcage coil (accuracy of 10% or better); (v) temperature rise vali- 
dation in a phantom with the actual long femoral nail implant (accuracy of 20% 
against the experiment) and; (v) temperature rise validation in the detailed human 
model with the same femoral nail implant located at approximately the same depth 
(accuracy of 25% or 2.4 °C on average with the standard deviation of or 0.2 °C 
against the experiments with the phantom). 

In the last case, the tool testbed predicted the higher temperature rise (by approx- 
imately 25% or 2.4 °C higher on average) at the implant tips than the in vitro experi- 
ments with the simplified gel phantom. An additional validation of the MDDT was 
therefore made against in vivo measurements in living human arm which indicated 
the temperature deviation of the MDDT from the in vivo experiment of only 10%. 

This tool addresses the growing number of patient MRI scans and expected prev- 
alence of patients with implanted medical orthopedic devices specifically mid-aged 
and elderly women who are most affected by osteoporosis and the associated bone 
fracture. Newly developed implants, as well as legacy medical implants without 
MRI safety information, need to be evaluated for safety in the MRI environment. 

The widely used ASTM F2182 standard allows measuring RF implant heating in 
a homogeneous gel phantom. However, the human body anatomy is different from 
this phantom. For safety and completeness, it may be therefore desired to addition- 
ally estimate the heating of the same implant in a realistic heterogeneous anatomical 
human model. 

This MDDT tool may help to estimate the induced RF heating of a metallic 
orthopedic implant in the detailed model of an elderly woman aged ~60 at the fre- 
quency of 64 MHz (1.5 T MRI scan). This tool can also help to identify the appro- 
priate worst-case implant size, configuration, orientation, and an allowable scan 
protocol (coil power) by performing multiple simulations to determine the 
RF-induced temperature rise as a function of the scan protocol and the required 
scan time. 

This tool has demonstrated that it accurately predicts absolute temperature rise 
for RF-induced heating with acknowledgement of the following limitations: 


* Although all separate blocks of the modeling pipeline (human model, coil model, 
SAR values) were validated separately and in the general case, the end result — 
temperature distribution along an implant and as a function of time — is limited 
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to one yet most critical implant geometry: a long femoral nail (a nearly straight 
metal rod) subject to the most excessive heating and possible resonant effects. 

* The focus of this tool is currently at 1.5 T/64 MHz cylindrical bore MRI systems 
using circularly polarized RF birdcage full body coils. 

* The tool currently does not include multi-parts orthopedic implants or implants 
with screws. 

* The tool currently does not provide heating estimates for gradient-coil induced 
power deposition for orthopedic implants. 

e The computational human model margins are limited to women aged 50-70 and 
with a high obesity (BMI of 30-36). 


The main advantage of using this tool lies in the possibility to accurately estimate 
temperature rise near an orthopedic implant in the realistic high-resolution elderly 
female human model and thus additionally justify and assess measurements obtained 
with the ASTM F2182 standard and their safety margins. 

While modeling SAR distributions near implants in realistic virtual human mod- 
els is well understood, accurate modeling and prediction of temperature rise at and 
near implants in human models is much more difficult. The reason is a necessity to 
couple an electromagnetic software and a thermal software with a realistic hetero- 
geneous human model including fine implant and tissue details which are best 
described with the finite element method that is developed for curved and fine 
geometries. 

The present tool addresses this knowledge gap by constructing and validating the 
state-of-the-art multiphysics FEM modeling pipeline. It couples the detailed accu- 
rate CAD human model, the world-best Ansys electromagnetic modeling software 
HFSS including the MRI coil modeling, and the Ansys thermal software, all within 
one single user-friendly shell — the Ansys Workbench. 

The second advantage of using this tool is the type of the embedded anatomical 
female CAD model that is appropriate for elderly female subjects of 50—70 years 
old with a high obesity. There is no similar anatomical CAD model currently 
available. 

The main disadvantage with the tool is that it currently requires long computa- 
tion times necessary to obtain accurate and reliable temperature rise results. These 
times are on the order of 6—24 hours depending on the computer hardware used. 

As stated previously, the tool's context of use is currently restricted to the long 
femoral nails and does not include other types of the femoral and shoulder orthope- 
dic implants. It is also is limited to the 1.5 T cylindrical bore MRI system only. 
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2 Tool Description and Validation 


2.1 VHP Female v 5.0 Anatomical Human Model (DOI 
10.20298/VHP-FEMALE-V.5.0) 


The VHP Female v.3.0—5.0 computational human model (also known under the 
nickname ‘Nelly’ to users of Dassault Systemes SIMULIA software CST), shown 
at left in Fig. 1, is an anatomically accurate heterogeneous female (~60 years old, 
~88 kg, BMI of ~36, obese, with heart pathology) surface-based human body model. 
It was constructed from the photographic cryosection data of the Visible Human 
Project of the US National Library of Medicine with the world-best isotropic resolu- 
tion of 0.33 mm. Its construction has been well documented in the literature [5—10], 
but several points are worth emphasizing here: 


1. There are 249 distinct components or triangular 2-manifold surface meshes 
(with an additional 40 characterizing embedded implants). No intersections or 
joint faces between discrete meshes are allowed, enabling unique assignment of 
electromagnetic, thermal, or other material properties. 

2. The source data for the model are freely and publicly available [11, 12]. The 
complete co-registration data for all model cross-sections are also made publicly 
available [13]. 


One advantage of the model is its topological and computational simplicity. The 
model size is purposely limited to approximately 0.4 M triangular facets in total 
with an attempt to keep the anatomical accuracy within the body within 1-7 mm. 
This makes it possible to apply virtually any commercial or custom-tailored finite- 
element or boundary element computational solver in a very reasonable amount of 
time. Nearly identical results are obtained when using different software packages 
and numerical methods [14]. Yet another advantage is the separation of the muscu- 
lar system of the body into individual muscle objects (about 50 in total). 

Some recent uses of the VHP Female v3.0—5.0 computational human model have 
been documented in the literature — cf. for example [14—31] — and compared with 
experimental data on radio frequency propagation. A number of studies have 
focused on the characteristics and behaviors of antennas near the human body [15- 
23]. Others have examined through simulation various biomedical applications 
including transmission channel modeling [24], transcranial direct current stimula- 
tion [25], estimation of bone density [26], gastroenterology [27, 28], SAR simula- 
tions in MRI coils at 3 T [29], and safety of active implants under MRI procedures 
[30, 31]. 

The VHP-Female model has been licensed by several private entities engaged in 
medical device modeling including Ansys, Inc., Dassault Systémes SIMULIA, 
Stryker Corp., Cambridge Consultants, Inc., and WIPL-D. A simplified version of 
the model, VHP-Female v. 2.2, has over 500 registered college users and is available 
for free download online. 


Computational Tool Comprising Visible Human Project® Based Anatomical Female... 137 


z= 80 mm grid :10 mm per division 

l WhiteMatter 
Skull 
Skin 
GrayMatter 
FatShell 
CSFShell 


Veinsupper 
AvgBodyShell Veins lower 
Uina Radius right 


Ulna. Radius eft 


Stomach 
Spleen 
——— SkinSiell 
Saas <= " Ribs.right12 
~ Ribs.right11 
« BB Ribs righti0 


Ribs right9) 
Ribs.right& 


Ribs.elt12 
Ribs leftt1 
Y Ribs left10 


Ribs left9 


Bg Rosiers 
Muscle Rectus Abdominal right.top 
Muscle Rectus Abdominal left top 
Muscle.LatissimusDorsi right 
Muscle.LatissimusDorsi left 
Muscle. Forearm Flexors.right 
Muscle Erector spinae right 

E 


Muscle Forearm Extensors.right 


Muscle. Forearm Flexors left 
Muscle Forearm.Extensors left 
Muscle Erector spinae left 
Muscle.Bisep right 

Muscle Bisep left 

Humerus right 
SpinalCordGrayMatter 
FatShell 

SpinalCordCSF- 
CartilageG.right 

Cartilages left 

 AvgBodyShell 

Arteries 


Veins.lower 

SkinShell 
Sciatic,Nerve.right 
Sciatic.Nerve.left 
Muscle. Sartorius.right 
Muscle. Sartorius.left 
Muscle.Quadriceps.right 
Muscle.Quadriceps.left 
Muscle.Hamstring.right 
Muscle.Hamstring.left 
Femur.right 

Femur.left 
Femur.Bone.Marrow.right 
Femur.Bone.Marrow.left 
FatShell 

AvgBodyShell 

ll Arteries 


Fig. 1 (a) Full-body Visible Human Project Female v5.0 CAD based computational phantom 
embedded into Ansys FEM software. The phantom is composed of 249 individual structures. Some 
individual muscles are removed for clarity. (b-e) Examples of co-registration maps for surface 
meshes in different transverse planes with tissue labeling. In case (b), non-anatomical separation 
between scalp and skull was corrected. In case (d), the tissue labeling list is not shown. The com- 
plete full-body co-registration maps with the vertical resolution of 1 mm and with tissue labeling 
in every cross-sectional plane are available online in *.mp4 format [13] for independent inspection 
and verification purposes 


2.2 Ansys FEM Computational Software Suite 
2.2.1 Computation of SAR 


The human model has been coupled with the proven Ansys Electronics Desktop 
(HFSS) finite-element electromagnetic software with adaptive mesh refinement. It 


138 G. Noetscher et al. 


has been employed to solve Maxwell’s equations in three-dimensional space. In this 
way, the model is considered as an arbitrary (inhomogeneous) isotropic medium 
with piecewise-constant electric permittivity e having the units of F/m and with 
constant magnetic permeability u having the units of H/m. 

After Maxwell's equations for an electric field (or the electric field intensity) 
E(r,t) [V/m] and for a magnetic field (or the magnetic field intensity) H(r, t) [A/m] 
are solved, the local SAR (W/kg) is defined through averaging the dissipated power 
per unit mass over a small (ideally infinitesimally small) volume V, that is 


sin) JE Qr av T 


where oc is (generally piece-wise constant) medium conductivity with the units of 
S/m, p(r) is the local mass density, and IE(r)l is the electric field magnitude at the 
observation point. The body-averaged or the whole-body SAR;,;, is given by averag- 
ing over the entire body volume, as 


1 o(r) 2 
ee E(rY av 5 
m 5. d Ipe) (r) 2) 


Similarly, SAR;, is given by averaging over a volume with the weight of 1 g 


SAR, (r) -- | Deef dV (3) 


SARjo,(r) is found in a similar fashion. 


2.2.2 Computation of Temperature Rise 


The SAR is the density of volumetric heat sources, not the temperature rise itself. 
To compute the resulting local temperature rise, the solutions of electromagnetic 
simulations have been coupled to Ansys Mechanical/Thermal via Ansys Workbench. 
Pennes' bioheat equation, based on the heat diffusion equation, is a standard approx- 
imation for heat transfer in biological tissues [39-41] implemented in the Ansys 
Mechanical software. It has the form 


p(r)C(r) E =V{k(r) V7) =p(r) 0+ p(r)san- 
p.CB(T-T,)=0 e 
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where T(r,t) is the local temperature, p(r) is density, C(r) is specific heat capacity 
(at constant pressure), Q is the metabolic heat generation rate, B is the perfusion rate 
[1/s]; index b is related to blood. 

All electromagnetic, mechanical, and thermal tissue properties used in the pres- 
ent study are catalogued in the reputable IT'IS database [42]. The electromagnetic 
material properties are given as a function of frequency, which is ideal for the pres- 
ent analysis. 


2.3 Validation of Overall Segmentation and Anatomical 
Correctness of the Human Model 


2.3.1 Topological Validation 


All 249 distinct components of the human model have been proven to be strictly 
2-manifold (watertight) surface meshes or a reasonable mesh quality. All intersec- 
tion or joint faces between discrete meshes have been eliminated enabling unique 
assignment of electromagnetic, thermal, or other material properties. 


2.3.2 FEM Software Compatibility Validation 


The entire model (including the implant parts) was so constructed that is passes the 
optional built-in Ansys mesh intersection and mesh quality checker with the highest 
option "Strict". Such a check can be independently performed and confirmed by any 
MDDT user. 


2.3.3 Anatomical Validation 


To validate the model anatomically, the co-registration method has been used, which 
implies direct superposition of transverse cross-sections of all surface tissue meshes 
onto the original cryosection images. An in-house MATLAB module (Fig. 1b-e) 
was written that performs such superposition with the resolution of 1 mm and 
simultaneously labels all tissue meshes which are present for a given cross-section. 
Its output is a scanning sequence in *.mp4 format [13]. 

After resolving multiple mesh intersections, the validation of the model segmen- 
tation shown in Fig. 1a was performed by 


(i) co-registration of surface meshes superimposed onto the original cryo-section 
images (Fig. 1b-e) and; 
(ii) correcting all intersection flaws. 


Further visual anatomical validation was performed by a number of anatomical 
experts in their respective areas from Beth Israel Deaconess Med. Ctr. and 
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Massachusetts General Hospital, Boston MA including Profs A. Nazarian (ortho- 
paedics), A. R. Opotowsky, (cardiovascular systems), V. Poylin, (gastroenterology), 
E. K. Rodriguez (nusculoskeletal tissue components), A. Pascual-Leone (cranial 
and intracranial anatomy) and Prof. G. Haleblian (urology). The questionable sur- 
face meshes (mostly bones but also soft tissues including scalp, bladder, uterus, etc.) 
were corrected. An example of a non-anatomical flaw is shown in Fig. 1b where the 
scalp was non-anatomically separated from the skull during the cryogenic process, 
which required proper mesh adjustment. 

The surface deviation for all meshes was found not to exceed 1—8 mm. The high- 
est accuracy was achieved for 


(i) extracerebral and intracranial volume; 
(ii) vertebral column/spinal cord. 


The complete full-body co-registration maps (movies) with the vertical resolution 
of 1 mm and with tissue labeling in every cross-sectional plane are available online 
to the MDDT users in *.mp4 format [13] for independent inspection, verification 
and design/development purposes. 


2.4 Validation of Overall SAR Prediction by the MDDT During 
a 1.5 T MRI Scan 


The first step of the MDDT workflow — the computation of local SAR — has been 
validated by comparison of normalized SAR predicted by two different modeling 
techniques in a 1.5 T birdcage whole body MRI coil. 

A generic, whole-body high-pass birdcage coil with 16 rungs and 32 matching 
capacitors, loaded with the VHP Female model, has been considered. The coil has a 
diameter of 64 cm and length of 69 cm consistent with [46]. The simulation geom- 
etry is shown in Fig. 2a. Simulations have been conducted at shoulder/heart and 
abdominal landmarks. For the former, the coil center is oriented to coincide with the 
top of the T7 vertebra; the latter has the coil center located at the top of the L1 ver- 
tebra. The coil was tuned to the desired frequency of 64 MHz (B, = 1.5 T) when 
loaded with the subject at each landmark (cf. [53]). Similar to [53], an ideal excita- 
tion was applied with 32 sources placed in the two end-rings to perform the function 
of the capacitors. This excitation provides results which are very similar to the con- 
ventional two-port or four-port excitations [53]. 

A comparison was further made with the results of Ref. [53] which was using a 
nearly identical high-pass birdcage coil (diameter of 63 cm and length of 70 cm), 
the nearly identical heart landmark, and the identical coil excitation type. However, 
an in-house voxel model for the same VHP-Female dataset was employed in [53] 
followed by the FDTD simulation method with the resolution of 5 mm. 

To compare the solution variation as a function of FEM mesh density and the 
solution convergence trends, solutions were generated first with 1 and then with 8 
adaptive mesh refinement passes, created approximately 0.5 M and 2.0 M tetrahedra, 
respectively. As an example, Fig. 2 shows the corresponding local SAR distributions 
at two different FEM resolutions in the coronal plane for the coil loaded with the 
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Fig. 2 Local SAR distribution in the coronal plane for a high pass full-body RF coil operating at 
64 MHz loaded with the VHP-Female v5.0 computational phantom given B, amplitude of 1 pT 


at the coil center. (a) Positioning of model within the birdcage at shoulder/heart landmark. (b) 
Ansys HFSS (Electronics Desktop) solution with one adaptive pass. (c) Ansys HFSS (Electronics 
Desktop) solution with eight adaptive passes 


VHP-Female v3.0 computational human model given B, amplitude of 1 uT at the 
coil center. Fields solutions from Ref. [53] were also normalized given the desired 
magnitude of B, at the coil center of 1 uT. The normalization is done in the form 


SAR 
SAR > ———— (6) 


(B; AuT) 


The local SAR was computed in Ansys HFSS and then exported to MATLAB over 
a uniform 3D grid of 2 mm in size. Whole-body SAR;,;, was computed from this 
data directly in MATLAB. SAR;, was also calculated by finding a volume surround- 
ing the observation point having the mass of exactly 1 g, and then performing aver- 
aging according to Eq. (3c). This averaging volume contains approximately 
5 x 5 x 5 individual voxels (2 x 2 x 2 mm each) closest to the observation point. The 
observation points form a 3D sub-grid spaced of 20 mm and 10 mm, respectively. 
The tissue density was set at 1 g/cm? uniformly in space. SAR,,, was computed in 
the same way. In this case, the averaging volume contained approximately 1250 
individual voxels (2 x 2 x 2 mm each) closest to the observation point. 

Figures 2b, c show the local SAR within the various tissues of the VHP-Female 
v3.0 model after 1 pass and 8 passes, respectively, for the shoulder/heart landmark. 
All results are normalized to the B, amplitude of 1 uT at the coil center. By ana- 
lyzing Fig. 2 we can see that 


1. The accurate solution with eight adaptive mesh refinement passes generates a 
more realistic SAR distribution, especially with regard to the local SAR — see 
Table 1. In particular, two non-physical maxima of the SAR observed at the top 
of the head are no longer present. 
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Table 1 Comparison of SAR values predicted using the VHP-Female surface CAD model with the 
values predicted by the reference voxel model derived from the identical image dataset [53]. 
Nearly identical high-pass birdcage coil dimensions, coil landmark, and excitation type were used. 
All results are normalized to 1 pT Bj field at the coil center 


Max. Max. Max. 
Whole- |non- lg 10g 
Coil body averaged local local 
Source | Method Model landmark SAR local SAR SAR SAR 
Present | FEM - 1 CAD Shoulder/ 0.16 44.5 5.22 2.61 
report | adaptive VHP-fem. | heart (top of 
pass (Ansys) | v. 3.0 vert. T7) 
88 kg 
64 MHz 
Present | FEM- 8 CAD Shoulder/ 0.13 12.0 1.61 1.37 
report | adaptive VHP-Fem. | heart (top of 
passes v. 3.0 vert. T7) 
(Ansys) 88 kg 
64 MHz 
Ref. FDTD VoxelVis. | Heart 0.12 NA 1.78 NA 
[53] Voxel size Human 
5mm Female 
64 MHz 


2. The accurate solution with eight adaptive mesh refinement passes and the coarse 
solution with one adaptive pass generate approximately the same whole-body 
SAR and SAR distribution maps, but considerably different peak local SAR val- 
ues — see Table 1. Note that the coil has been retuned separately in both cases. 

3. The maximum SAR for the present landmark is observed in the upper shoulder/ 
neck area and in the arms area. 

4. A gap between the arms and the body may generate large SAR values. 


Table 1 compares the computed SAR values with the values obtained in Ref. [53] 
under nearly identical conditions. For both models and both methods, the whole 
body SAR and the maximum SAR,,(r) differ by 8% and 10%, respectively. 

A comparison with a variety of other modeling reference sources [35—37, 46—56] 
has also been made both at 1.5 T and 3 T. It was found that the SAR results calculated 
when using the VHP Female v3.0 model are strictly within the bounds of all other 
reported values. 


2.5 Validation of Temperature Rise Prediction by the MDDT 
Near a Long Femoral Nail Implant Against Measurements 
with an ASTM-Like Phantom 


ASTM and ASTM-like phantoms have been used as the standard method for testing 
compliance of implants within MRI environment [54]. These phantoms are not an 
equivalent of intricate human body structures. Here, the comparison is made with 
the original experiments performed in a series of papers [43-45]. 
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2.5.1 Modeling Testbed 


Prior to conducting simulations using the VHP Female v5.0 model, a necessary 
workflow and RF power calibration have been established. Figure 3 depicts to scale 
simulation of the original experiment performed in a series of papers [32-34]. It was 
also used for initial calibration purposes. At left in Fig. 3, a homogenous experimen- 
tal AGAR gel phantom with electromagnetic properties matching that of the human 
body is positioned within a 1.5 T birdcage MRI coil operating at 64 MHz. 

This phantom (40 cm long by 20 cm wide) in Fig. 3 is not the exactly the ASTM 
F2182 phantom. However, it does support and even enhances the loop-like distribu- 
tion of the electric field in the birdcage coil and thus describes a realistic heating 
scenario, perhaps even the worst-case scenario with respect to implant heating. 

All examinations [32-34] were performed with a 1.5 T MR scanner 
(MAGNETOM Symphony, SIEMENS) with the phantom at the center of the coil. 
A birdcage shaped transmit/receive body coil was used there with the inner diame- 
ter of 60 cm and the length of 70 cm, which is close to the present modeling setup 
(64 cm and 69 cm, respectively). Coil tuning and excitation was performed as 
described in the previous example. We again use ideal excitation which is very simi- 
lar [53] to the two-port excitation used in [32-34]. 


Simulated Current Density 


Simulated Temperature Rise 


Fig. 3 Simulations to establish calibration for experiments conducted in [32-34]. At left, the 
metallic nail implant placed within a homogeneous loop-like Agar phantom [32-34] at a depth of 
2 cm. At top right, the current density produced with a whole-phantom SAR of 4.0 W/kg. At bot- 
tom right, the simulated temperature given a total volumetric power loss of 120 W exactly corre- 
sponding to experiment [32-34]. Simulation results produced a temperature rise of 10.02 ?C, 
slightly less than the 12.6 °C experimentally observed in [32-34] 
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A 24 cm long metallic orthopaedic nail implant of Zimmer, Inc. made of non- 
magnetic stainless still (originally very slightly bent but modeled as a straight rod of 
the same diameter) has been used. Such implants are normally used to treat frac- 
tures in clinical practice. 

The implant was embedded in this phantom such that it is 2 cm away from the 
top and side edges of the phantom. Both the phantom and metallic rod have been 
assigned dimensions and material properties consistent with experiment [32-34]. 


2.5.2 Replication of Original Experiments 


The computed volumetric current density within the phantom is shown at top right 
in Fig. 3. This density matches well with the published result; it was produced by 
adjusting coil power to exactly replicate the RF exposure given in [32-34] — a 
whole-phantom SAR of 4.0 W/kg and a volumetric power loss density of 120 W in 
the phantom. The temperature simulation is shown at bottom right of Fig. 3. Given 
the above power loss density as an internal heat generation source, the total tem- 
perature rise within the human model was 10.02 ?C. This value is slightly less than 
the 12.6 °C observed in [32-34]. However, it is close enough (a 20% temperature 
difference) to give confidence in the adequate simulation setup, enabling the exten- 
sion of this methodology to the case involving the VHP female human model. 


2.5.3 Heating of the Same Implant in the MDDT Human 
Female Phantom 


Figure 4 at left shows the VHP Female computational human model in the 1.5 T 
birdcage MRI coil. The human model has been oriented so that the center of the 
femur bone is aligned with the center of the coil. The left quadriceps muscle within 
the VHP model is not shown in Fig. 4 so that the position of the femur can clearly 
be seen. Within the left femur, the same 24 cm long cylindrical metallic implant has 
been inserted as required in the surgical practice. 

This rod is assigned material properties consistent with [32-34]. The coil is driven 
at 64 MHz and at a power such that the volumetric power loss within the human 
model is again 120 W, also consistent with the published experimental results. 

A length-based mesh constraint of no edge larger than 2.5 mm was enforced for 
the metallic rod mesh and a total of about 636,000 tetrahedral elements were used 
in the Ansys HFSS simulation. 

Once the electromagnetic simulation was complete, the results were passed to 
Ansys Mechanical software by linking the two simulations in Ansys Workbench. A 
mesh refinement was again employed on the faces of the rod to ensure that a dense 
enough mesh was created to capture the local temperature changes. Approximately 
629,000 tetrahedral elements were used in the transient thermal simulation. Volumetric 
losses (SAR) produced by the Ansys HFSS simulation were imported into Ansys 
Mechanical and used as the internal heat generation source density. The RF coil and the 
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Fig. 4 The VHP-Female computational phantom positioned with the 1.5 T MRI birdcage coil. 
Some body parts are removed for clarity. At left, the femur position is illustrated to show its orien- 
tation within the model — the metallic nail implant is aligned to reside within the trabecular bone 
structure. At top right, each individual object within Ansys Mechanical model is assigned specific 
thermal properties. At bottom right, the temperature rise is shown after 900 seconds of continuous 
coil operation. These values correspond well with published experimental data [32-34] 


heat sources were active for 900 s and the model was allowed to cool for another 600 s. 
The values for the implant temperature captured from 0 to 900 s are shown in Fig. 5. 

Figure 5 shows temperature dynamics: it compares simulated numerical values 
(red stars) obtained using the VHP Female computational phantom with published 
experimental data (two black curves) [32-34] for the maximum temperature rise 
near the implant. The depth of the implant within the VHP model is approximately 
3-4 cm and the maximum temperature value after 900 seconds of coil operation 
time is represented by the top red star. 

Peak simulated temperature values for the implant within the human model are 
shown in the bottom right of Fig. 4. A maximum temperature rise of approximately 
11 °C is computed at the very ends of the metallic implant against 8.25 °C observed 
in experiment at 4 cm depth [32-34]. Thus, the MDDT testbed predicts a higher (by 
2.75 °C) maximum temperature rise of 8.25 °C than what was measured in [32-34] 
which constitutes the difference of 3396 or 2.75 ?C. 

This is likely due to the different (non-homogeneous) material properties 
employed in the present study, the realistic angled orientation of the nail, and a pos- 
sible resonant behavior of the long implant. In the VHP Female model, the rod is 
aligned with the femur and represents a more realistic position that would be 
encountered in a clinical setting. 
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Fig. 5 Comparison of the simulated numerical values (red stars) obtained using the VHP Female 
computational phantom with published experimental data (two black curves) [32-34] for the maxi- 
mum temperature rise near the implant. The depth of the implant within the VHP model is approxi- 
mately 4 cm and the maximum temperature value after 900 seconds of coil operation time is 
represented by the top red star. The VHP model predicts a slightly higher (by 2.75°) maximum 
temperature rise than what was measured in [32-34]. This is likely due to the different (non- 
homogeneous) material properties employed in the present study and the slightly angled orienta- 
tion of the nail 


2.6 Variation and Uncertainty of Temperature Rise 
Measurements Near the Implant 


Table 2 gives the least mean square difference in the temperature rise (maximum 
implant temperature at the tip is recorded) over the entire interval from 0 to 15 min 
between the modeled data and the in-vitro data at 2 and 4 cm implant depth, respec- 
tively, in degrees C. Statistical significance of the deviation is quantified via the 
p-value for the paired-sample t-test with p < 0.05 considered statistically 
significant. 


2.7 Validation of Temperature Rise Prediction by the MDDT 
Against In Vivo Measurements 


In the previous case validation example, the MDDT testbed predicted the higher 
maximum temperature rise (up to 33% higher) at the implant tips than the experi- 
ment in vitro with the simplified gel phantom. An additional validation of tempera- 
ture rise was therefore made against in vivo measurements in living human. 

Since in-vivo measurements in the MRI coil at 1.5 T have not been reported, a 
comparison was made with tissue heating due to a small single-loop coil close to 
skin surface at 165 MHz [38]. The replication of the experimental setup [38] is 
shown in Fig. 6. The forearm of the VHP Female computational human model has 
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Table 2 Least mean square difference in the temperature rise (maximum implant temperature at 
the tip) over the entire interval from 0 to 15 min between the modeled data and the in-vitro data at 
2 and 4 cm depth in degrees C 


Deviation: 2 cm Deviation: 4 cm Deviation: 2 cm Deviation: 4 cm 
Phantom implant depth vs. implant depth vs. |implantdepth vs. | implant depth vs. 
vs. MDDT  |leftleg left leg right leg right leg 
23^C 2.5 °C 23°C 256 
Side view Bottom view 
a) 


Fig. 6 The forearm of the VHP-Female computational human model adjacent to a single loop coil 
operating at 165 MHz in Ansys HFSS. The model geometry is shown, including internal forearm 
tissues, each defined with specific electromagnetic and thermal properties 


been isolated to speed up the simulation. Internal structures, including the humerus, 
ulna and radius bones, extensor, flexor, triceps and biceps muscles, radial nerve, and 
various arteries and veins, are encapsulated within concentric layers of muscle, fat 
and skin tissue. The coil [38] is modeled as a 80 mm diameter copper torus with a 
minor diameter of 2 mm driven at 165 MHz by a 50 Ohm lumped antenna port and 
a thin sheet of Teflon separates the antenna from the forearm. These settings exactly 
correspond to experiment. All thermal and electromagnetic material properties 
associated with the internal body structures correspond to those supplied in [38]. 

According to [38], the coil antenna is provided sufficient power to dissipate 
approximately 31 W within the forearm and produce peak local SAR values in the 
lower corners of the forearm of approximately 450 W/kg. These power conditions 
have been replicated in simulations with the forearm of the VHP Female model. 
Adaptive mesh refinement produced a volumetric mesh consisting of approximately 
0.1 M tetrahedra. 

Following simulation in Ansys HFSS, the peak SAR value experienced in the 
outermost model layer is computed as approximately 450 W/kg, consistent with the 
value reported in [38] to within 5%. 

Then, the geometry and resulting fields were passed via Ansys Workbench to 
Ansys Mechanical. Here, the human model was again remeshed, producing just 
over 0.06 M elements. The volumetric power losses determined in Ansys HFSS 
were imported and evaluated as the heat generation source density. Heat generation 
sources were active for 120 s and the model then cooled via convection for 19.2 s; 
temperatures were recorded throughout the duration of the simulation. The corre- 
sponding temperature rise after 139 seconds in total was computed and compared 
with experiment. 
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The simulated and measured [38] thermal maps within the forearm in a plane 
directly above the coil center have been examined. The peak simulated temperature 
change is about 6.6 °C, which is slightly higher that the measured 6 °C change 
reported in [38]. The deviation is within 10%. 

The spatial temperature map are also consistent with the measured map despite 
the different wrist geometry and tissue composition. There is some difference in the 
observed depth of temperature change. This is likely mostly due to the fat layer of 
the VHP Female model which is thicker than the subject reported in [38]. 


3 Conclusions 


Based on the evidence provided, this non-clinical assessment tool entitled 
“Computational Tool Comprising Visible Human Project® Based Anatomical 
Female CAD Model and Ansys HFSS/Mechanical® FEM Software for Temperature 
Rise Prediction near an Orthopedic Femoral Nail Implant during a 1.5 T MRI Scan” 
was found to reliably predict both temperature distribution and its evolution in time 
along the long femoral nail metal implant caused by RF power deposition from the 
1.5 T birdcage MRI full body coil. The tool can also help to identify the appropriate 
worst-case device and coil size, configuration, and orientation by performing mul- 
tiple simulations to determine the RF-induced temperature rise as a function of a 
scan protocol and required scan time. 

All separate blocks of the modeling pipeline — the human model topology and 
anatomy including co-registration, surgically correct implant embedding, the RF 
coil model, and the resulting SAR and temperature behavior — have been validated 
independently and all the validation results have been made available to the user. 
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High Frequency Electromagnetic 
Modeling: Microwave Imaging 


Modeling and Experimental Results A 
for Microwave Imaging of a Hip gese 
with Emphasis on the Femoral Neck 


Johnathan Adams, Peter Serano, and Ara Nazarian 


1 Introduction 


Osteoporosis affects approximately 21.2% of women and 6.3% of men over the age 
of 50 world-wide [1]. In the United States alone, the estimated economic burden of 
osteoporosis-related fractures in 2005 was $17 billion and is expected to increase to 
about $25 billion by 2025 [7]. Hip fracture is one of the most serious and debilitat- 
ing outcomes of osteoporosis [2, 3], with a 14-36% mortality rate during the first 
year following a fracture [4]. Hip fracture incidence rates increase exponentially 
with age in both women and men [5]. In 2010, there were estimated to be 158 mil- 
lion persons at high risk for a bone fracture, a staggering statistic expected to double 
by 2040 [6]. 

The World Health Organization (WHO) has defined individuals at risk for these 
fractures based on their areal Bone Mineral Density (aBMD, g/cm’) relative to that 
of a normal young adult, as measured by Dual-energy X-ray Absorptiometry (DXA) 
[8]. However, DXA is not without its flaws. These include exposing patients to 
small ionizing radiation doses of up to 0.86 mrem [9]; measurement errors due to 
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surrounding soft tissues [10, 11]; bone mineral density (BMD) measurements are 
affected by variations in bone size [12, 13]; and cortical and trabecular bone cannot 
be differentiated [14]. Additionally, fracture predictions based on aBMD are neither 
sensitive nor specific [15—19]. While DXA-based aBMD has been previously shown 
to be an important predictor of hip fracture risk [20], it does not offer a direct assess- 
ment of bone’s load bearing capacity [21, 22]. Additionally, the predictive BMD 
value for fracture decreases in individuals over 70 years [23]. From the age of 60 to 
80, the risk of hip fracture increases 13 times, while the decrease in BMD can only 
account for doubling of the risk [24]. Also, there is a wide overlap in BMD scores 
of postmenopausal women [25] who do and do not sustain osteoporotic fractures 
[26], and approximately 50% of fragility fractures occur in patients with DXA- 
derived BMD T-scores in the normal or low bone mass range [27, 28]. 

Quantitative ultrasound devices provide a low-cost, non-ionizing alternative to 
DXA using a specialized ultrasound transceiver to measure bones near the surface 
of the skin. A commercial example of this, Bindex®, uses the pulse-echo technique 
to measure the thickness of the frontal cortical shell of the tibia bone [29-32]. It has 
been found to correlate strongly with DX A measurements [29], a less than optimal 
gold standard. 


1.1 Why Microwave Imaging 


Microwave or radiofrequency imaging of (heel) bone was first introduced by Dr. 
Keith Paulsen and his research group at Dartmouth College in 2010, as an alterna- 
tive non-ionizing diagnostic method to assess bone health [2, 33—36]. Due to the 
well-known complexity and poor spatial resolution of the standard microwave 
imaging setup [37, 38] used in these studies, they produced no clinically applicable 
results. However, the underlying physical idea of this method is simple and power- 
ful. We have previously designed a simpler device to prove the viability of the con- 
cept at the wrist [39] and have achieved approximately 83% sensitivity and 94% 
specificity using a neural network classifier to differentiate between osteoporotic 
and healthy subjects [40]. 


1.2 Potential Difficulties 


Taking these radiofrequency measurements is not without challenges. The transmis- 
sion must pass through the bone in the region of interest and arrive at the receiver 
antenna with sufficient power to be measured, given a range of individuals. This 
makes measuring bones deeper in the body more difficult compared to more super- 
ficial bones. Additionally, the antennas must be placed such that the major compo- 
nent of the received signal is through the bone rather than its surrounding tissues. 
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1.3 Our Approach 


This study consists of a set of simulations to determine field propagation inside the 
body validated by in vivo experimental measurements under the same conditions. 
The simulations produced models that included reflection coefficient 5;, and trans- 
mission coefficient $5, in addition to the fields. These S-parameters can be measured 
in a physical setup using a network analyzer. The simulations and physical measure- 
ments were performed with the same antennas [41]. Additional simulations were 
performed with different antennas to investigate wideband measurements; these 
were not verified experimentally. The simulation results were analyzed primarily 
based on the electric field and Poynting vector. 


2 Materials and Methods 


This study was divided into two parts: first, a set of in vivo measurements using real 
antennas and second, a set of simulations using a corresponding human body model. 
The measurements were taken with Institutional Review Board (IRB) approval 
(IRB-19-0123) through Worcester Polytechnic Institute. The same human subject 
was used for all in vivo measurements. 


2.1 Experimental Hardware 


The antennas featured in this study are dual antiphase patch antennas [41] built 
using copper on FR4. Two sets of antennas, shown in Fig. 1, were investigated. 

Set A (resonators: 2.0 cm x 1.4 cm, ground-plane: 5.0 cm x 1.9 cm) connected 
to matching networks that match them to 675 MHz. Matching networks were built 
with lumped components and applied at the antenna feeds, after the 180° power 
splitter (Mini-Circuits® ZFSCJ-2-232-S+, 5 MHz to 2.3 GHz). 

Set B (resonators: 2.5 cm x 1.6 cm, ground-plane: 3.0 cm x 8.0 cm) were not 
matched to any particular frequency. The antenna feeds connected directly to the 
180° power splitter (Mini-Circuits® ZFRSC-183-S+, DC to 1.8 GHz). 

Both antennas had 0.5 cm spacing between the resonators. The antennas were 
connected to a Keysight FieldFox N9914A network analyzer. The network analyzer 
transmitted at —15 dBm over a frequency range of 30 kHz to 2 GHz at 401 points. 
The magnitude in dB and phase in degrees of S|; and $5, were saved to a CSV-file. 
The measurements were each a single frequency sweep. 
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Fig. 1 Comparison of size between antennas from Set B (top) and from Set A (bottom). The left 
two are the physical antennas and the right two are the corresponding CAD models. The spacing 
between patches in both antennas is 0.5 cm. Antennas in Set A had resonators of 2 cm x 1.4 cm and 
a ground-plane of 5 cm x 1.9 cm. Antennas in Set B had resonators of 2.5 cm x 1.6 cm and a 
ground-plane of 3 cm x 8 cm. The antennas were fed from the back, the solder joints in the figure 
are the feeds 


2.2 Measurement Sites 


To test the viability of various sites for measuring transmission through the femoral 
neck, we first checked using both sets of antennas to determine if meaningful trans- 
mission could occur given the positions of the antennas. The exact positions inves- 
tigated are shown in Fig. 2. 

The positions investigated were: 


. On the side of the body, positioned over the greater trochanter. 

. On the side of the body, positioned next the iliac crest. The antenna in this posi- 
tion was rotated in the plane of the drawing in Fig. 2 to investigate different 
polarizations. The orientation shown in the figure (vertically aligned with the 
body and the antenna in position 1) was considered 0^, and rotation angles were 
measured toward the front of the body (clockwise on the right side, counter- 
clockwise on the left side). 

. On the front of the body, positioned over the anterior superior iliac spine. 

. On the rear of the body, positioned over the top edge of the gluteus maximus. 

5. On the front of the body, positioned horizontally in the same horizontal plane as 

the greater trochanter. 

6. On the rear of the body, positioned horizontally and below the gluteus maximus. 


Ne 


RW 


Measurements were taken between two of the positions. The antennas were held to 
the body by the subject being measured, by pressing on the center of the ground 
plane of each antenna. This ensured deformation of the body so that the total length 


Modeling and Experimental Results for Microwave Imaging of a Hip with Emphasis... 159 


A) Coronal Plane B) Sagittal Plane 
2 4 
pelvis 3 
1 (/ 
5 
femur / 


Fig. 2 (a) Front (coronal plane) view of the right side of the body with antenna positions by num- 
ber. (b) Right side (sagittal plane) view of the body with antenna positions noted by number. In 
both, the skin profile, right pelvis, and right femur are shown in addition to the antennas. Antennas 
were pushed against the body such that gaps, such as the one near position 4, were not present 
during the measurement. The antenna in position 4 was located over the posterior iliac crest 


of the antenna was contacting skin. The positions that were not in use for a given 
measurement did not have antennas present. All position combinations measured 
were measured with both Set A and Set B antennas. 


2.3 Simulated Antenna Positioning and Human Body Model 


Antenna positions on the simulated body model were the same as those on the in 
vivo model and are shown in Fig. 3. The base CAD model is the Ansys male human 
model. It was chosen to match the in vivo subject, who is male. 

The CAD model includes the full body with bones, muscles, and fat modelled 
throughout. Some skin layers and fat deposits are represented in aggregate by a 
volume with the average electric properties of the human body [42, 43]. Cartilage in 
joints, such as the hip joint, is not modelled by default. We investigated the effects 
of cartilage by producing a new shell using the area between the femur and pelvis 
making up the ball joint. This new volume was between 3 to 10 mm thick due 
mostly to the large-triangle tessellation of the bones' shells. In addition, two outer 
skin shells were added with properties derived from the VHP-Female v.5.0 
model [44]. 
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Fig. 3 (a) Front (coronal plane) view of the right side of the CAD model with antenna positions 
by number. (b) Right side (sagittal plane) view of the CAD model with antenna positions noted by 
number. In both, the wireframe body shell, right pelvis, and right femur are shown in addition to 
the antennas. The body shell was flattened or Boolean-subtracted using the antenna's shape to 
eliminate gaps and ensure good coupling at each position. The apparent difference in location of 
position 4 is due to perspective of the drawing in Fig. 2 


2.4 Software Modelling of Matched Antennas 


The matching networks were modelled in Ansys using S-Parameter measurements 
of the physical matching networks. The matching networks’ measurements were 
taken over the same frequency range (30 kHz to 2 GHz) and with the same resolu- 
tion (401 points) as the in vivo transmission measurements. However, they were 
taken at 0 dBm and averaged across 8 sweeps, whereas the in vivo measurements 
were taken at a lower power (—15 dBm, see above) and with only a single frequency 
sweep. The input ports were all 50 Q characteristic impedance, identical to the 
physical network analyzer. The 180? power splitters were modelled with an ideal 
splitter model. Figure 4 shows a typical simulation configuration for a single antenna 
at position 1 compared to a physical measurement at the same site. 
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Comparison Between Simulated and Measured Reflection Coefficients 
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Fig. 4 Comparison of reflection coefficient magnitude IS, ,| for the simulated antenna (red) and the 
in-vivo antenna (dashed black). This figure additionally shows the configuration of the matching 
and power splitting circuits in HFSS for the simulated curve. Both the simulated and measured 
curves were produced from antennas at position 1 


3 Results 


This section is divided between the in vivo and simulated results. The in vivo results 
provide an assessment of the total transmission in a real subject, given a set of 
antenna positions, and the simulations show where inside the body the transmis- 
sion occurs. 


3.1 In Vivo Measurements 


The positions shown in Fig. 2 are positions between which transmission was 
achieved. Additional sites were measured, including one between sites 4 and 6 
through the center of the gluteus maximus, but no meaningful signal was received. 
Figure 5 shows the transmission coefficient for a selection of antenna position pairs 
using Set A, while Fig. 6 shows the transmission coefficient for the same position 
pairs using Set B. 

In addition to the measurements shown in the figures, transmission from position 
1 to position 2 was measured with varying polarizations achieved by rotating the 
position 2 antenna in 45? increments. For Set A, the highest average transmission 
over the bandwidth of the antenna transmission was seen at 45? of rotation while the 
lowest at 135°. Set B was less consistent but showed similar results: minimum 
transmission at 270? and maximum at either 45? (on the left side of the body) or 
135? (on the right side). Overall, Set B showed less change in $5, over various 
angles of rotation than Set A, potentially due to noise. Set A showed differences of 
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Set A, S5 for Various Paths 
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Fig. 5 Comparison of transmission coefficient $5; when using antennas from Set A for three prop- 
agation paths: First, position 1 to position 2, a semicircular path through the compartment. Second, 
position 5 to position 6, a straight path through the upper femur. Third, position 3 to position 4, a 
straight path through the upper pelvis 


$5; at the same frequency within the passband of the antenna with different polariza- 
tions up to 20 dB, while Set B showed differences up to only 10 dB. Lying down 
during the measurement process decreased this variation by about half. Table 1 
shows the maximum measured magnitude of the transmission coefficient for each 
orientation tested, for both sets of antennas. 


3.2 Simulations 


First, the relative agreement between the simulated and measured results is charac- 
terized by Fig. 7, in which there are resonances at approximately the same frequen- 
cies in the measured and simulated environments, but the simulated environment 
experiences significantly more attenuation on transmission than the measured 
environment. 

Next, the propagation paths of the waves were observed using animated electric 
field plots in various observation planes and 3-D Poynting vector plots in the bones 
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Set B, S5 for Various Paths 
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Fig. 6 Comparison of transmission coefficient $5, when using antennas from Set B for three prop- 
agation paths: First, position 1 to position 2, a semicircular path through the compartment. Second, 
position 5 to position 6, a straight path through the upper femur. Third, position 3 to position 4, a 
straight path through the upper pelvis 


Table 1 Maximum transmission coefficient magnitude |S,,| for various relative polarizations 
using both sets of antennas on the right leg. Measured from position 1 to position 2 


Set A, right leg Set B, right leg 
Angle, degrees f. MHz 155,1, dB f. MHz 155,1, dB 
0 (co-polarized) 615.0 —38.663 10.0 —12.254 
45 635.0 —29.127 10.0 —12.066 
90 (cross-polarized) 695.0 —42.658 10.0 —13.860 
135 705.0 —39.581 10.0 —11.955 
180 (co-polarized) 705.0 —37.817 10.0 —12.741 
225 600.0 —36.861 10.0 —12.891 
2770 (cross-polarized) 710.0 —39.151 10.0 —14.000 
315 705.0 —37.694 10.0 —12.057 


and the body. At higher frequencies, a surface-propagating wave is present, as seen 
in Fig. 8. The Poynting vector plots in the femur and pelvis for the three transmis- 
sion configurations in Figs. 5, 6, and 7 are shown in Fig. 9. 
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Comparison of Transmission Coefficient for 3 Paths 
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Fig.7 Comparison of transmission coefficient $5, when using antennas from Set A for three prop- 
agation paths: First (red): position 1 to position 2, a semicircular path through the compartment. 
Second (blue): position 5 to position 6, a straight path through the upper femur. Third (green): 
position 3 to position 4, a straight path through the upper pelvis. The dashed lines are the measured 
in-vivo IS,,| (also seen in Fig. 5) and the solid lines are simulated 


Fig. 8 Electric field magnitude in the sagittal plane at different frequency bands. (a) is 60 MHz, 
(b) is 550 MHz, and C is 715 MHz. All three are snapshots from animations, taken at a phase of 
60°. Note the vertical surface-propagating wave is present in b and c but not in a. The antennas for 
these measurements are the Set A antennas, located at position 1 


In addition to the results shown in the figures, simulations were performed with 
a dielectric “belt” between the transmitting and receiving antennas to attenuate the 
surface wave. The effect was not strong enough to reduce the magnitude of the sur- 
face wave to a level comparable to that of the wave propagating through the bone. 


4 Discussion 


4.1 Limitations 


This study only considered the strongest component of the received wave. 
Simulations suggest that this component likely propagates through skin, fat, and 
muscle when the antennas are on the same side of the body compartment. The same 
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Fig. 9 Poynting vector distribution in the femur and pelvis for three antenna position pairs: (a) 
transmission from position 1 to position 2, (b) transmission from position 3 to position 4, and (c) 
transmission from position 4 to position 5. Poynting vector magnitude is represented by color, 
warmer is larger 


simulations also suggest propagation occurs through the bone and this second com- 
ponent will accrue some phase shift (delay) relative to the one that propagates 
through the soft tissue. 

This study performed in vivo measurements on only one subject, a 26-year-old 
male, who is not at significant risk of osteoporosis according to standard risk factors. 

In vivo spectra were determined from a single frequency sweep; therefore, noise 
contents in the spectra are more significant than had the measurements been per- 
formed using averaging of multiple sweeps. 

Matched (Set A) and unmatched (Set B) antennas are not identical and have dif- 
ferent resonant frequencies. Set A had a bandwidth of about 230 MHz, centered on 
about 675 MHz (when matched) and Set B had a bandwidth of about 420 MHz, 
centered on about 215 MHz. 


4.2 Validation of Simulation Using In Vivo Results 


While the measured and simulated results are not a perfect match, the resonant fre- 
quencies are collocated in the two spectra for the same antenna set and positions. 
Some of the difference in transmission coefficient magnitude between the measured 
and simulated spectra is due to differences between the model and the physical 
subject. These differences include the level of detail of the CAD model, and differ- 
ences in physical shape between the CAD model and the subject. 
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4.3 Propagation Paths and Antenna Position 


To achieve transmission through the bone, antennas should be placed on the oppo- 
site sides of the body compartment. If placed on the same side of the body compart- 
ment, the surface wave has the shortest path between the two antennas and thereby 
dominates the received signal. Direct transmission across a body compartment, con- 
trarily, puts the shortest path between the antennas through the bone at the center of 
the compartment, and the Poynting vector for such a setup is the largest at the center 
of the compartment [39]. This is illustrated in Fig. 9, where the distribution in part 
C shows more even transmission through the femoral neck than part A. Part C’s 
antennas are transmitting across the compartment while part A’s transmit in a 
u-shape, starting and ending on the same side of the compartment. 


4.4 Frequency Choice for Propagation Through Bone 


Itis common knowledge that lower frequencies provide better human body penetra- 
tion but lower spatial resolution in microwave imaging. Our simulations have con- 
firmed this, but we also note that in this application we can consider frequencies that 
are lower than traditionally considered for microwave imaging of the human body, 
due to the independence of this approach from spatial features. Therefore, a fre- 
quency of operation closer to 60 MHz with a reduced surface wave is preferable in 
this application to a higher frequency. Any waveguide-like effects from higher- 
frequency waves propagating through bones are overshadowed by the lack of pen- 
etration to reach the bones in the first place, and by the large surface-propagating 
waves produced by these high frequencies. 
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